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ABSTRACT 


The  purpose  of  this  work  is  to  establish  how  a  moving  emitter  can  be  localized  by  a  passive 
receiver  through  the  use  of  out-of-plane  multipath  signals  reflected  by  the  terrain.  This  is  a  novel 
localization  technique  that  assumes  no  a  pviori  knowledge  of  the  location  of  the  multipath  sources.  The 
emitter  parameters  of  range,  heading,  velocity,  and  altitude  are  estimated  by  exploiting  the  correlation 
between  the  direct-path  signal  and  the  delayed,  attenuated,  and  Doppler  modulated  signals  reflected  by 
the  terrain. 

Two  basic  assumptions  about  the  terrain  scattering  properties  lead  to  different  maximum 
likelihood  (ML)  estimators  of  emitter  parameters.  The  first  assumption  is  that  the  terrain  scattering  is 
fundamentally  homogeneous;  in  this  case  the  ML  estimator  is  found  to  have  the  structure  of  a  time- 
varying  FIR  filter.  The  second  assumption  is  that  the  terrain  scattering  is  fundamentally  inhomogeneous 
and  dominated  by  a  number  of  discrete  scatterers.  This  assumption  leads  to  a  two-part  estimator  which 
first  estimates  the  scattering  parameters  of  azimuth,  differential  delay,  and  Doppler  shift,  and  then 
estimates  the  emitter  parameters  using  the  scatterer  parameter  estimates.  The  Cramer-Rao  lower  bounds 
of  the  variances  for  each  estimator  are  derived  and  used  to  study  estimator  performance  for  several 
scenarios. 

The  proposed  estimators  are  successfully  demonstrated  using  field  data  collected  at  White  Sands 
Missile  Range  during  the  DARP A/Navy  Mountaintop  program.  Several  extensions  to  the  basic  results  are 
considered  such  as  localizing  pulsed  and  self-correlated  emitters,  multiple  emitters,  and  the  effects  of 
receiver  motion. 
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The  material  presented  in  this  report  is  identical  to  that  in  a  dissertation  submitted  to  Northeastern 
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1.  INTRODUCTION 


1.1  BACKGROUND 

Passive  localization  of  an  emitting  source  is  a  problem  that  has  received  considerable  interest  from 
both  the  radar  and  sonar  communities.  Localization  refers  to  the  estimation  of  source  parameters  such  as 
the  range,  angle,  heading,  velocity,  and  altitude  or  depth  of  a  radiating  source  by  receiving  and 
processing  the  source  emissions.  For  military  applications,  a  passive  localization  system  offers 
advantages  over  an  active  one  because  it  operates  covertly.  That  is,  it  can  perform  the  localization 
without  active  transmissions  that  could  give  away  its  own  position.  Strictly  speaking,  a  radar  operating  in 
a  passive  mode  is  not  acting  as  a  radar  but  as  a  radio  frequency  (RF)  receiver.  Hence,  when  the  term 
radar  is  used  in  this  report,  we  are  usually  referring  to  the  radar  receiver. 

Classes  of  emitters  that  will  be  considered  include  communications  transmitters,  pulse 
transmitters,  and  noise  transmitters  or  jammers.  A  jammer  is  a  high  power  emitter  that  radiates  a  noise¬ 
like  signal  designed  to  disable  a  radar  or  sonar  system  by  increasing  the  noise  level  [1].  Very-low 
sidelobe  or  adaptive  antennas  can  be  employed  to  reduce  the  susceptibility  to  jamming  via  the  antenna 
sidelobes  [2]-[5];  however,  if  the  jammer  is  in  the  active  sensor’s  mainbeam,  the  sensor  may  be 
completely  disabled.  Thus,  an  active  sensor  is  not  able  to  determine  the  range  to  a  jammer  because  it 
cannot  detect  its  own  relatively  weak  echoes  traveling  the  same  propagation  path  as  the  very  strong 
jamming  signal.  The  jammed  sensor  will  be  able  to  determine  the  jammer  azimuth  and  elevation  angles 
by  passively  receiving  and  performing  angle-of-arrival  estimation  (AOA),  or  direction  finding  (DF), 
using  the  jammer’s  strong  direct-path  signal.  However,  these  angle  measurements  alone  do  not  yield  the 
range  to  the  jammer.  Since  emitters  of  this  type  often  represent  important  targets,  techniques  to  localize 
these  emitters  are  important  to  both  the  radar  and  sonar  communities. 

The  work  presented  in  this  report  is  primarily  geared  towards  radar  with  the  Mountaintop  UHF 
phased-array  radar  as  the  principal  sensor.  Extensive  Mountaintop  data  containing  airborne  emitters  have 
been  collected  for  the  analysis  of  terrain  scattered  interference  (TSI)  [6], [7]  (or  so-called  “hot  clutter 
[8]”)  and  bistatic  clutter  [6], [9], [10].  This  same  data  will  be  used  to  study  emitter  localization  techniques 
presented  herein  [1 1],[12].  The  potential  use  of  TSI  for  emitter  localization  was  suggested  first  by  Gerard 
Titi  [11]. 

Before  introducing  the  novel  out-of-plane  multipath  localization  concept  that  is  the  subject  of  this 
report,  we  first  review  some  basic  multisensor  and  single-sensor  localization  techniques.  These  are 
presented  to  introduce  the  localization  problem  and  to  define  some  commonly  used  localization  terms. 

1.1.1  Multiple  Sensor  Localization  Techniques 

There  are  numerous  passive  techniques  that  can  be  used  to  localize  an  emitting  source.  For 
radio-wave  emitters,  [13]  briefly  discusses  nineteen  passive  ranging  techniques,  and  for  underwater 
acoustic  sources,  [14]  discusses  several  ranging  and  direction  finding  techniques.  A  few  of  these 
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techniques  are  discussed  in  this  section  in  order  to  introduce  some  localization  terms  and  concepts.  They 
provide  background  for  the  single-sensor  techniques  described  in  the  next  section  and  the  out-of-plane 
multipath  technique  introduced  in  Section  1.2. 

One  of  the  simplest  forms  of  localization  is  triangulation  which  is  shown  in  Figure  1.1  (a)  and  (b). 
In  part  (a),  two  physically  separated  but  cooperating  receivers  obtain  angle-of-arrival  estimates  (AOA), 
which  are  also  called  direction  lines,  or  lines-of-position  (LOP),  by  measuring  the  angle  to  the  emitter 
using  some  angle  measurement  technique.  (A  discussion  of  AOA  and  DF  techniques  can  be  found  in 
[15]-[17].)  The  master  station  receives  angle  measurements  from  the  other  stations  in  the  localization 
network,  over  data  links,  and  determines  the  emitter  position.  Cooperating  AWACS  (ESA)  radars  utilize 
a  passive  mode  to  localize  emitters  in  this  manner  [18].  If  more  than  one  emitter  is  present,  false 
associations  may  occur  and  generate  false  targets  that  are  commonly  referred  to  as  ghosts.  This  is 
illustrated  in  Figure  1.1  (b)  for  a  two  station,  two  emitter  scenario.  Ghost  targets  can  often  be  removed  by 
employing  cross-correlation  and  time-difference  of  arrival  (TDOA)  techniques  [19],[20]. 


Station  1  Station  2 

(Master) 


(a) 


(b) 


Figure  1.1.  (a)  Emitter  localization  via  triangulation,  (b)  Two  emitters  generate  four  intersections  resulting  in  two 
ghost  or  false  target  positions  in  addition  to  the  true  target  positions. 
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Hyperbolic  measurements  systems  use  TDOA  measurements  instead  of  angle  measurements  to 
localize  an  emitter.  An  example  using  three  isotropic  sensors  and  an  emitter  in  a  horizontal  two- 
dimensional  plane  is  shown  in  Figure  1.2.  The  master  station  receives  replicas  of  the  emitter  waveform 
from  the  other  stations  over  a  data  link.  TDOA  estimates  are  obtained  by  cross-correlating  the  emitter 
signals  received  at  the  various  stations.  The  time  differences  specify  hyperbolic  curves  and  the 
intersection  of  these  curves  yields  the  emitter  position.  For  localization  in  three  dimensions,  four  sensors 
are  commonly  used  and  the  emitter  position  is  taken  from  the  intersection  of  three  hyperboloids.  If  there 
is  relative  motion  between  the  emitter  and  the  sensors,  differential  Doppler  measurements  may  also  be 
utilized.  Hybrid  localization  systems  combine  multiple  discriminants,  such  as  TDOA,  AOA,  or 
differential  Doppler. 
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Figure  1.2.  Three  sensor  hyperbolic  localization  scheme. 
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1.1.2  Single  Sensor  Localization  Techniques 

Single-sensor  localization  techniques  rely  on  virtual  sensors  to  perform  localization  estimates.  The 
virtual  sensors  are  formed  by  multiple  propagation  paths  (multipath)  or  by  sensor  motion.  Once  the 
virtual  sensors  are  realized,  localization  can  be  performed  using  multisensor  techniques  like  those 
described  in  the  preceding  section.  The  obvious  advantage  to  the  single-sensor  techniques  is  that  only  a 
single  sensor  is  required;  no  multisensor  coordination  or  data  links  are  necessary.  The  single  passive 
sensor  functions  independently.  Passive  localization  using  out-of-plane  multipath,  the  topic  of  this  report, 
is  a  single  sensor  technique. 

Sensor  motion,  such  as  that  achieved  with  an  airborne  radar  receiver,  permits  AOA  measurements 
to  be  performed  at  multiple  sensor  positions.  The  AOA  measurements  are  used  to  triangulate  the  emitter 
in  the  same  manner  as  the  multisensor  technique  illustrated  in  Figure  1.1(a).  Since  the  sensor 
measurements  are  performed  sequentially,  this  technique  is  often  called  sequential  localization.  This 
method  works  well  when  the  receiver  is  moving  and  the  emitter  is  stationary  but  may  have  difficulties 
when  the  emitter  is  moving.  In  contrast,  we  shall  see  in  Chapter  5  that  the  performance  of  the  estimators 
designed  in  this  report  improves  with  increasing  emitter  velocity. 

Several  in-plane  multipath  techniques  may  be  used  to  perform  single-sensor  localization.  A  few 
simple  examples  are  shown  in  Figure  1.3.  This  figure  shows  a  radar  using  reflections  from  the  ground  or 
ionosphere,  and  a  sonar  using  reflections  from  the  ocean  floor  or  surface,  to  localize  the  emitter.  The 
virtual  sensors  are  the  images  of  the  true  sensor  generated  by  multipath.  Autocorrelation  of  the  sensor 
signal  yields  TDOA  estimates  between  the  direct  propagation  path  and  the  reflected  paths.  Angle 
measurement  techniques  can  be  used  to  measure  the  angle  of  arrival  of  the  direct  and  reflected 
propagation  paths.  The  TDOA  or  AOA  measurements,  when  combined  with  knowledge  of  the  position  of 
the  reflecting  surfaces,  are  used  to  perform  the  localization  using  the  multisensor  techniques  described  in 
the  preceding  section. 

Knowledge  of  the  position  of  reflecting  surface  is  essential  for  these  single-sensor  localization 
techniques.  That  is,  the  high-frequency  (HF),  over-the-horizon  (OTH)  radar  receiver  needs  to  know  the 
height  of  the  ionosphere,  and  the  general  radar  receiver  using  ground  reflection  needs  to  know  its  altitude 
relative  to  the  reflection  point.  Similarly,  the  sonar  needs  to  know  its  depth  and  distance  to  the  ocean 
floor.  Errors  in  knowledge  of  the  position  of  the  reflecting  surfaces  translate  directly  to  errors  in  emitter 
localization.  The  out-of-plane  multipath  technique  developed  in  this  report  requires  no  prior  knowledge 
of  scatterer  location. 
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Figure  1.3.  Single  sensor  localization  using  multipath  (assumes  a  flat  earth  model). 


1.1.3  Background  Literature 

The  preceding  two  sections  provide  a  brief  qualitative  introduction  to  localization  using  multiple- 
and  single-sensor  techniques.  The  reader  is  referred  to  the  very  rich  literature  on  emitter  localization  if 
more  information  is  required.  Some  of  the  papers  we  found  helpful  are  indicated  in  this  section. 

Most  of  work  in  the  emitter  localization  literature  discusses  one  or  more  of  the  four  basic 
components  of  the  problem:  the  geometrical  model,  the  statistical  model,  parameter  estimation,  and  the 
accuracy  or  performance  assessment.  The  geometrical  model,  including  sensor  configurations,  is  unique 
for  a  specific  problem.  The  equations  describing  the  emitter  and  source  geometries  for  localization  are 
usually  nonlinear,  such  as  the  intersection  of  hyperbolas  discussed  in  Section  1.1.1.  The  statistical  model 
usually  describes  measured  data  collected  by  the  sensors  in  the  presence  of  measurement  noise.  For 
example,  when  measurement  noise  is  present,  more  than  two  hyperboloids  will  not  intersect  at  a  single 
point.  The  problem  becomes  that  of  detennining  the  most  likely  position  and  assessing  the  accuracy  of 
the  result.  The  parameter  estimation  component  is  concerned  with  determining  the  best  emitter 
parameter  set  based  on  some  optimization  criterion,  such  as  least-square  error  or  a  maximization  of  the 
likelihood  function.  The  accuracy  or  performance  assessment  usually  contains  a  detemunation  of  the 
Cramer-Rao  lower  bounds  (CRLB)  which  specify  the  minimum  variance  an  unbiased  estimator  can 
achieve.  If  sub-optimum  parameter  estimators  are  employed,  they  are  usually  accompanied  by  an 
evaluation  of  bias  and  a  comparison  to  the  CRLB. 

A  discussion  of  statistical  models  and  the  other  components  of  the  localization  problem  can  be 
found  in  [21]-[23]  for  many  of  the  multisensor  localization  techniques.  These  references  linearize  the 
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nonlinear  geometrical  equations  by  expanding  in  a  Taylor  series  and  retaining  only  the  linear  terms. 
These  results  are  expanded  to  include  sensor  position  uncertainty  in  [24], [25].  Single  sensor  localization 
techniques  using  reflection  from  the  ionosphere  are  described  in  [26],[16],[27]-[29].  These  techniques 
can  be  thought  of  as  vertical  triangulation  [26],  and  the  localization  sensors  are  called  single-site  locators 
(SSL).  Sequential  localization  techniques  for  stationary  emitters  and  moving  sensors  are  found  in 
[20], [30].  Sequential  techniques  that  also  employ  Doppler  shift  are  found  in  [31],[32].  The  accuracy  of 
systems  using  virtual  sensors  generated  by  multipath  is  discussed  in  [33]  for  two  sensors  and  two  paths 
(direct  and  reflected)  and  in  [34]  for  an  arbitrary  number  of  sensors  and  paths.  The  sonar  practice  of 
matched  field  processing  (MFP)  [35]-[37],  is  a  technique  whereby  a  plane-wave  beamformer  is  replaced 
by  a  generalized  beamformer  with  weighting  vectors  computed  from  a  multipath  model  of  the  ocean 
waveguide.  These  techniques  localize  a  source  in  range,  depth,  and  azimuth,  and  may  include  as  many  as 
forty  or  more  rays  or  modes  in  the  multipath  model.  It  is  considered  an  inverse  problem  that  is  solved  by 
constructing  an  “infinite”  number  of  forward  solutions  on  a  parameter  grid  and  selecting  the  best  match 
to  the  observed  signal. 

Localization  methods  that  utilize  relative  time  delay  or  Doppler  frequency  difference  estimates 
(i.e.,  TDOA  or  FDOA  estimates)  can  draw  upon  the  extensive  literature  on  these  subjects.  The  maximum 
likelihood  estimator  for  time-delay  estimation  is  based  on  cross-correlation  and  is  described  in  [38],[39]. 
Reference  [39]  is  a  tutorial  review  of  coherence  and  time-delay  estimation  which  contains  an  additional 
74  references.  Time-delay  estimation  using  LMS  adaptive  tapped-delay-line  filters  is  described  in 
[40],[41].  Estimation  of  differential  Doppler  shift  is  discussed  in  [42].  The  maximum  likelihood 
estimator  for  the  joint  estimation  of  differential  delay  and  Doppler  is  discussed  in  [43]  where  it  is  shown 
that  the  estimation  of  differential  Doppler  is  uncoupled  from  differential  delay  and  phase  estimators.  A 
discussion  of  localization  using  both  TDOA  and  FDOA  methods  for  a  stationary  transmitter  and  two 
moving  receivers  can  be  found  in  [44], [45]. 

1.2  INTRODUCTION  TO  LOCALIZATION  USING  OUT-OF-PLANE  MULTIPATH 

This  section  introduces  the  concept  of  localization  using  out-of-plane  multipath  that  is  the  subject 
of  this  report.  This  is  a  new  localization  technique  that  makes  no  a  priori  assumptions  about  the  locations 
of  the  scatterers.  The  acoustic  (sonar)  and  electromagnetic  (radar)  multipath  examples  depicted  in  Figure 

1.3  were  examples  of  in-plane  multipath  localization.  The  direct  propagation  path  and  the  forward 
scattered  path  are  contained  in  a  vertical  plane  perpendicular  to  the  horizontal  reflecting  surfaces.  These 
are  specular  reflections  from  flat  surfaces  that  obey  Snell’s  law  of  reflection,  i.e.,  the  angle  of  incidence 
equals  the  angle  of  reflection.  Out-of-plane  scattering  refers  to  reflections  from  scatterers  in  the 
horizontal  plane  and  is  illustrated  in  Figure  1.4,  which  shows  a  receiver,  emitter,  and  reflections  from 
three  scatterers.  These  reflections  are  much  weaker  than  the  in-plane  reflections  described  above  but  for  a 
strong  emitter  such  as  a  jammer,  they  may  be  appreciable.  The  fact  that  these  reflections  are  significant  is 
evidenced  by  the  interest  in  the  terrain  scattered  interference  (TSI)  problem  [6],[8]-[10].  These  papers 
describe  the  use  of  adaptive  signal  processing  techniques  to  remove  the  TSI  signals  from  the  desired 
radar  signal.  The  unwanted  jammer  reflections  described  in  those  papers  are  the  same  signals  that  will  be 
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exploited  here  to  perform  emitter  localization.  This  is  yet  another  instance  of  the  adage  “One  radar’s 
clutter  is  another  radar’s  target  [46].’’ 


Figure  1.4.  Top-down  two-dimensional  view  of  the  radar,  emitter,  and  three  scatterers  [s,,  S2,  S3]  in  the  horizontal 
plane.  (A  three-dimensional  view  of  this  geometry  is  provided  in  Figure  2.2.) 


Before  developing  the  out-of-plane  multipath  localization  concept  further,  we  need  to  make  some 
assumptions  about  the  capabilities  of  the  receiver.  The  receiver  must  be  capable  of  performing  direction¬ 
finding  or  angle  estimation  in  the  horizontal  plane.  That  is,  the  receiver  is  capable  of  measuring  the 
angle-of-arrival  of  the  emitter  and  azimuthal  angle,  of  the  scatterers.  The  angle  of  the  emitter  is  used 
to  orient  the  coordinate  system  and  is  thus  set  to  zero  degrees  in  this  section.  In  addition  to  estimating  the 
angle  of  the  scatterers,  the  receiver  must  also  be  capable  of  estimating  the  time  difference  of  arrival 
(TDOA),  Xi,  between  the  emitter  direct  path  and  the  scattered  paths.  As  we  shall  see  below,  the  receiver 
must  also  be  capable  of  measuring  the  Doppler  frequency  difference  of  arrival  (FDOA),/,  for  the  direct 
and  scattered  paths.  The  Mountaintop  system  [6]  meets  these  requirements.  This  system  has  a  nominal 
bandwidth  of  200  kHz,  a  beamwidth  of  8  degrees,  and  a  center  frequency  of  435  MHz. 


To  proceed  with  our  development,  we  first  consider  the  case  of  TDOA  measurements,  [Ti,  T2,  x{\, 
and  angle  measurements,  [^1,  ^],  for  the  three  scatterers  and  a  stationary  emitter  shown  in  Figure  1.4. 

Can  we  localize  the  emitter  using  this  information?  Each  scatterer  forms  a  bistatic  triangle  with  sides  a„ 
bi,  and  R,  where  R  is  the  desired  localization  quantity  of  emitter  range.  The  TDOA  is  T;  =  (a,  +  fo,  -  R)/c, 
where  c  is  the  speed  of  light.  The  TDOA  estimates  are  converted  to  range-difference  measurements  by 
multiplying  by  the  speed  of  light,  Dj  =  (a,-  bi  —  R).  This  relation  and  the  law  of  cosines, 
bf  =  af  +  R^  -  la^R  cos  ({> I ,  can  be  solved  for  R  in  terms  of  the  range-difference  estimates  D„  and  range 
to  the  scatterer  a,-,  giving 


R  = 


°,-P,-P,V2 

D;  -  a.  -t-a,.  cos^,. 


(1.1) 
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Similar  equations  arise  in  bistatic  radar  problems  where  R  is  the  range  to  a  target  and  a,  is  the 
distance  from  the  bistatic  receiver  to  the  bistatic  transmitter  [47].  The  three  scatterers  give  rise  to  three 
equations;  however,  there  are  four  unknowns:  the  three  scatterer  ranges,  [a,,  02,  aj,  and  the  desired 
quantity,  R.  For  N  scatterers,  there  will  be  A^+1  unknowns  and  the  system  remains  underdetermined  with 
infinitely  many  solutions.  Therefore,  we  cannot  determine  the  range  to  this  emitter  without  a  priori 
knowledge  of  the  scatterer  locations.  This  is  demonstrated  graphically  in  Figure  1.5  which  shows  two 
bistatic  triangles  with  the  same  azimuth  and  range  difference  D  =  a^+b^—  =  Oj  +  ^2  ~  ^  • 

However,  if  the  range  to  a  single  scatterer,  a,  were  known,  along  with  the  measured  scatterer  azimuth  0 
and  range  difference  D,  the  emitter  range  could  be  determined  using  a  single  application  of  (1.1).  This  is 
the  case  where  a  priori  knowledge  of  the  scattering  points  is  available  and  is  similar  to  the  single-sensor 
techniques  described  in  Section  1.1.2.  This  case  will  not  be  explored  in  this  report. 


Figure  1.5.  A  determination  of  range  difference,  Di  =  (Oi  +  b,  -  R),  and  scatterer  azimuth  (j)  is  not  enough  to 
determine  the  range  to  an  emitter.  This  figure  shows  two  bistatic  triangles  with  the  same  range  difference,  D,  and 
scatterer  azimuth,  (j),  but  different  scatterer  ranges,  a„  result  in  different  emitter  ranges,  /?,. 


We  now  consider  the  case  of  an  emitter  moving  with  velocity  V  and  direction  or  heading  H.  The 
emitter  motion  produces  Doppler  shifts  in  the  scattered  and  direct-path  signals  that  provide  additional 
information  that  can  be  used  to  localize  the  emitter.  Specifically,  the  Doppler  shifts  provide  information 
that  may  be  used  to  determine  the  angles  [yi,  72,  73]  which  are  shown  in  Figure  1.6.  With  y;  determined, 
all  the  angles  of  the  bistatic  triangle  are  known  since  a.  -ISO  —  tj).  -7, and  was  measured  by  the 
receiver.  With  the  angles  known,  the  range  R  can  be  determined  using  the  law  of  sines  and  the  measured 
range  difference  D;  by  writing 


^siny.  ^  sin^; 
sin  a,,  sin  a. 
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Figure  1.6.  Top-down  view  of  the  radar,  emitter,  and  three  scatterers  [S],  S2,  S3].  Emitter  is  moving  with  velocity  V 
in  the  direction  H. 


The  receiver  cannot  measure  the  absolute  Doppler  shift  of  the  emitter  direct-path  signal  and  the 
individual  scattered  signals  because  it  does  not  know  the  original  emitter  waveform.  Instead,  the  relative 
Doppler  shift  between  the  direct  and  scattered  paths  is  obtained  using  cross-correlation  techniques.  The 
observable  quantity  is  the  Doppler  frequency  difference  (FDOA)  for  the  direct  and  scattered  signal  paths. 
Denoting  the  Doppler  shift  for  the  direct  path  as  fop  and  the  scattered  path  as/j,  the  differential  Doppler 
shift  isf=fs  -fdp-  With  measurements  of  differential  Doppler/  for  each  scatterer,  the  desired  angle  fi  can 
be  written  as 


7.  =  -H+cos  ^ 


(cos(H)- 


(1.3) 


where  A  is  the  wavelength  corresponding  to  the  center  of  the  receiver  narrowband  filter,  V  is  the  emitter 
velocity,  and  H  is  the  emitter  heading.  The  derivation  of  this  equation  is  contained  in  Section  2.4.2.  It  is 
derived  based  on  a  two-dimensional  approximation  that  ignores  emitter  altitude.  This  is  a  common 
approximation  for  bistatic  radar  problems  because  the  emitter  range  is  usually  much  greater  than  the 
emitter  altitude  [48].  A  more  general  three-dimensional  treatment  will  be  used  for  the  estimation  of 
emitter  altitude,  however. 


We  proceed  by  using  equation  (1.3)  in  (1.2)  to  yield 

^  A _ 

rsin7,(g,V,/.)  sin^.  ^ 

sina,.(i/,y,7.)  ^ 


(1.4) 
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which  expresses  the  emitter  range  R  in  terms  of  its  velocity  V,  heading  H,  and  the  measured  quantities  for 
each  scatterer  The  tilde  was  added  to  indicate  that  these  quantities  are  measured.  These 

equations  are  highly  nonlinear  and  must  be  solved  numerically.  Since  there  are  three  unknowns  [/?,  H,  V] 
at  least  three  scatterers  with  corresponding  estimates  of  differential  Doppler/,  differential  range  D„  and 
azimuth  0,  are  required.  In  practice,  a  large  number  of  scatterers  will  be  used  to  improve  estimation 
accuracy  and  to  eliminate  ambiguities  or  ghosts. 

In  this  section,  the  localization  using  out-of-plane  multipath  concept  was  introduced.  The 
technique  uses  a  combination  of  scatterer  AOA,  TDOA,  and  FDOA  measurements  resulting  in  a  set  of 
nonlinear  equations  that  can  be  solved  numerically.  This  approach  is  based  on  an  inhomogeneous  clutter 
assumption  because  it  relies  on  the  presence  of  a  number  of  dominant  scattering  points  and  distinct 
propagation  paths.  No  prior  assumptions  about  the  positions  of  the  scattering  elements  are  assumed, 
however.  A  homogeneous  approach  that  makes  different  assumptions  about  the  characteristics  of  the 
bistatic  clutter  will  be  developed  as  well.  The  basic  elements  of  the  localization  problem,  including  the 
geometrical  model,  statistical  model,  parameter  estimators,  and  an  accuracy  assessment,  will  be 
developed  for  both  the  homogeneous  and  inhomogeneous  clutter  models. 

1.2.1  Contributions 

The  fundamental  contribution  is  the  development  of  single-sensor  emitter  localization  techniques 
that  utilize  out-of-plane  multipath  signals  to  localize  moving  emitters  in  range,  heading,  velocity,  and 
altitude.  More  specific  contributions  include: 

1.  The  development  of  the  necessary  geometrical  relations  required  to  localize  emitters  using 
multipath. 

2.  The  development  of  statistical  models  for  the  description  of  emitter  multipath  signals. 

3.  The  design  of  estimators  for  the  localization  of  moving  emitters  using  out-of-plane  multipath  with 
no  a  priori  knowledge  of  scatterer  locations.  Both  homogeneous  and  inhomogeneous  scattering 
models  are  considered. 

4.  The  determination  of  estimation  accuracy  via  the  Cramer-Rao  lower  bounds  on  emitter  parameter 
variance. 

5.  The  demonstration  of  the  effectiveness  of  the  estimators  using  Mountaintop  field  data  containing 
airborne  emitters. 

6.  The  examination  of  several  extensions  to  the  basic  results,  such  as  the  localization  of  multiple 
emitters,  effects  of  receiver  motion,  and  performance  for  various  emitter  waveforms. 
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1.2.2  Report  Organization 

This  report  contains  eight  chapters  and  four  appendices.  Two  estimators  will  be  designed  that  are 
based  on  two  bistatic  clutter  modeling  approaches.  The  first  approach  is  to  assume  that  the  bistatic  clutter 
is  essentially  homogeneous  and  this  leads  to  one  estimator  design.  The  second  approach  is  to  assume  that 
the  bistatic  clutter  is  inhomogeneous  and  dominated  by  some  number  of  discrete  scatterers.  This  leads  to 
an  alternative  estimator  design.  Separate  chapters  will  be  devoted  to  the  development  of  each  estimator. 
The  contents  of  each  chapter  are  considered  below. 

Chapter  1  (this  chapter)  provides  background  on  general  localization  techniques  and  introduces  the 
out-of-plane  localization  concepts.  Chapter  2  extends  this  introduction  further  describing  in  more  detail 
the  bistatic  scattering  geometry  and  its  use  for  emitter  localization.  A  bistatic  Hough  transform  is 
developed  that  can  detect  an  intensity  locus  or  ridge  in  bistatic  delay-Doppler  space.  This  will  be  used  to 
examine  ambiguities  that  may  occur  using  this  localization  technique. 

In  Chapter  3  an  estimator  is  designed  based  on  the  assumption  that  the  bistatic  clutter  is  essentially 
homogeneous.  A  statistical  formulation  is  presented  and  the  maximum  likelihood  estimator  is  determined 
for  an  additive  Gaussian  noise  assumption.  The  estimator  is  found  to  have  the  structure  of  a  time-varying 
finite  impulse  response  (FIR)  filter.  The  complex  filter  coefficients  are  determined  using  least  squares. 
These  coeffieients  are  treated  as  nuisance  parameters  that  are  estimated  along  with  the  desired  emitter 
parameters.  Optimum  and  sub-optimum  estimation  techniques  will  be  presented. 

In  Chapter  4  an  estimator  is  designed  based  on  the  assumption  that  bistatic  clutter  is 
inhomogeneous  and  dominated  by  a  number  of  discrete  scatterers.  This  model  requires  a  two-part 
estimation  procedure.  The  receiver  first  estimates  the  pertinent  scattering  parameters,  such  as  relative 
time  delay,  r,  differential  Doppler,  /;  azimuth,  and  amplitude,  a,  for  some  number  of  the  dominant 
scatterers.  The  scatterer  estimates  are  then  used  to  estimate  the  emitter  parameters.  The  highly  nonlinear 
localization  equations  are  linearized  by  using  a  Taylor  series  expansion  and  retaining  only  the  linear 
terms. 

Chapter  5  contains  a  discussion  of  the  performance  of  the  estimators  designed  in  Chapters  3  and 
4.  The  Cramer-Rao  lower  bounds  (CRLB)  are  determined  and  evaluated  for  various  localization 
scenarios.  The  effect  of  varying  the  emitter  parameters  of  range,  heading,  velocity,  altitude,  and  power, 
as  well  as  the  receiver  parameters  of  bandwidth,  beamwidth,  steering  direction,  integration  time,  etc.,  will 
be  considered.  Monte-Carlo  simulations  are  performed  to  validate  the  bounds  calculations  and  overall 
implementation  of  the  estimators  prior  to  their  use  with  field  data. 

In  Chapter  6  we  use  the  estimators  designed  in  Chapters  3  and  4  to  localize  airborne  emitters 
contained  in  Mountaintop  data  collected  at  the  White  Sands  Missile  Range  (WSMR).  In  this  chapter  we 
localize  several  emitters  using  both  the  homogeneous  and  inhomogeneous  techniques.  Estimation  errors 
obtained  with  real  data  are  compared  to  the  CRLB.  This  chapter  serves  to  demonstrate  the  results 
presented  in  the  three  previous  chapters  using  experimental  field  data.  Problems  encountered  when 
applying  theoretical  results  to  experimental  data  will  be  discussed. 
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In  Chapter  7  we  show  how  the  basic  results  could  be  extended  to  include  a  broader  range  of 
scenarios  and  class  of  emitter  signals.  We  discuss  how  to  incorporate  Doppler  modulation  generated  by  a 
moving  receiver  as  well  as  that  generated  by  the  moving  emitter.  We  also  consider  the  localization  of 
multiple  emitters.  In  addition,  various  emitter  waveforms  are  examined  including  pulsed,  periodic,  and 
self-correlated  waveforms.  Generally  these  waveforms  cause  problems  for  localization  schemes  because 
multiple  peaks  arise  in  auto-  and  cross-correlation  functions.  The  multiple  peaks  may  generate  a  large 
number  of  false  emitters  or  ghosts.  However,  we  shall  see  that  the  out-of-plane  multipath  technique 
offers  some  immunity  against  these  problems  because  it  uses  FDOA  as  well  as  TDOA  information  in  the 
multipath  signals.  The  FDOA  information  arises  from  the  Doppler  shifts  generated  by  the  true  emitter 
which  may  be  effective  in  eliminating  ghosts. 

Chapter  8  contains  a  summary  and  conclusions  as  well  as  suggestions  for  future  research. 

Four  appendices  are  included  that  contain  the  necessary  mathematical  details  to  support  the 
estimator  designs  and  the  CRLB  calculations. 
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2.  LOCALIZATION  USING  OUT-OF-PLANE  MULTIPATH 


2.1  INTRODUCTION 

This  chapter  contains  additional  background  and  introductory  information  necessary  to  develop  the 
out-of-plane  multipath  estimation  techniques  that  are  the  subject  of  this  report.  It  begins  with  the  basic 
signal  and  noise  models  for  the  receiver.  Properties  of  the  emitter  waveform  and  receiver  noise  are  also 
described.  A  description  of  the  bistatic  scattering  geometry  and  its  use  in  solving  the  emitter  localization 
problem  is  also  provided.  The  chapter  concludes  with  the  development  of  a  bistatic  Hough  transform 
which  will  be  used  to  perform  emitter  localization  graphically  and  for  examining  potential  ambiguities 
that  may  arise  when  using  out-of -plane  multipath  for  emitter  localization. 

2.2  RECEIVER  SIGNAL  AND  NOISE  MODELS 

The  three  major  components  to  this  problem  are  the  emitter,  the  receiver,  and  the  ground.  These 
components  are  shown  in  Figure  2.1.  The  receiver  samples  the  emitter  waveform  directly  using  a 
reference  channel  and  the  scattered  replicas  of  the  waveform  using  the  multipath  channel.  One  or  more 
multipath  channels  may  be  used.  Each  channel  consists  of  a  directive  antenna  beam,  a  mixer,  a  lowpass 
filter,  and  an  analog-to-digital  (A/D)  converter.  The  directive  antenna  beams  permit  angle  measurement 
in  the  horizontal  plane.  No  directivity  in  the  vertieal  or  elevation  plane  is  assumed.  This  assumption  is 
consistent  with  airborne  surveillance  radars,  such  as  the  E2C  or  ESA,  which  generally  have  small  vertical 
apertures. 

The  received  signals  for  each  channel  are  amplified,  mixed,  filtered,  sampled,  and  stored  for 
processing.  The  amplifier  is  noisy  and  contributes  additive,  zero-mean,  white  Gaussian  noise  which  is 
independent  from  sample  to  sample  as  well  as  from  channel  to  channel.  A  quadrature  mixer  shifts  the 
center  frequency  from/o  to  zero-Hz  baseband.  A  lowpass  filter  band  limits  the  white  Gaussian  noise.  The 
bandwidth  of  the  lowpass  filter  is  matched  to  the  bandwidth  of  the  emitter  waveform  to  maximize  the 
signal-to-noise  ratio  of  the  received  signals.  The  result  is  a  narrowband,  eomplex  amplitude  or  complex 
envelope  representation  of  the  original  signal  received  at  the  antenna  terminals.  This  signal  is  sampled  by 
the  A/D  converter  above  the  Nyquist  rate  to  avoid  aliasing.  This  process  of  reducing  a  narrowband  signal 
to  its  complex  amplitude  representation  is  described  in  many  texts;  see,  for  example,  [4], [49]. 
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Emitter 


Receiver 


From  Ground 


Figure  2.1.  Basic  model  for  the  emitter,  receiver,  and  multipath.  The  reference  channel  collects  the  emitter  signal 
via  the  direct  path  and  the  multipath  channel  collects  the  signals  reflected  by  the  ground. 


2.2.1  Mountaintop  Radar  Receiver  and  Emitter 

The  localization  techniques  developed  in  this  report  will  be  demonstrated  with  field  data  collected 
as  part  of  the  DARPA/NAVY  Mountaintop  Program  [6].  This  program  utilized  the  14-element 
phased-array  radar  system  RSTER  (Radar  Surveillance  Technology  Experimental  Radar)  [50]  and  was 
designed  to  support  the  evaluation  of  space-time  adaptive  processing  (STAP)  algorithms  [51]  in  the 
presence  of  clutter  and  januning.  Jamming  was  provided  by  an  airborne  emitter  containing  a  2  kW 
transmitter  carried  by  a  Lear  Jet.  The  emitter  waveform  was  a  600  kHz  pseudo-noise  (PN)  sequence.  The 
phased-array  antenna  permitted  the  simultaneous  collection  of  multiple  receive  beams.  The  major  system 
parameters  are  listed  in  Table  2.1. 
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TABLE  2.1 


Mountaintop  Parameters  for  TSl  and  Emitter  Localization  Experiments 


Parameter 

Value 

Antenna  Gain 

29  dBi 

Antenna  Type 

14  El.  Phased  Array 

Antenna  Element  Gain 

17.5  dBi 

Operating  Frequency 

400-500  MHz 

Instantaneous  Bandwidth 

200  kHz 

Noise  Figure 

8dB 

Beamwidth  -  Azimuth 

8  degrees 

Beamwidth  -  Elevation 

6  degrees 

Radar  PRF 

312.5  Hz 

Max.  Coherent  Time 

153.5  ms  (48  PRls) 

Receiver  Altitude 

2.43  km  MSL,  -1 .2  km  AGL 

Airborne  Emitter  Power 

2  kW  over  600  kHz 

Emitter  Altitude 

~9.2  km  MSL,  -8  km  AGL 

2.3  BISTATIC  GEOMETRY  AND  SCATTERING  CONSIDERATIONS 

The  geometry  that  describes  bistatic  scattering  is  illustrated  in  Figure  2.2.  Three  angles  are 
required  to  describe  this  geometry.  The  angle  of  incidence  is  denoted  as  0i  and  the  scattered  angle  of 
incidence  is  denoted  as  6,.  The  out-of-plane  angle  is  denoted  as  and  is  defined  as  shown.  Instead  of 
incidence  angles,  the  complementary  depression  angles  are  sometimes  used.  Also,  some  investigators  use 
the  bistatic  angle  instead  of  the  out-of-plane  angle  when  describing  bistatic  scattering.  This  is  the  angle 
subtended  by  the  incident  and  scattered  rays. 
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Figure  2.2.  Bistatic  scattering  regions  and  geometry. 


^  scattering  regions  evident  in  Figure  2.2:  the  specular  point,  the  glistening  surface, 

and  the  diffuse  scattenng  region.  The  strongest  scattering  occurs  at  the  specular  point  (or  region)  which 
IS  represented  by  the  darkest  shading  in  the  figure.  At  this  point  the  out-of-plane  angle  is  zero  and  the 
angle  o^ncidence  is  equal  to  the  angle  of  reflection.  This  area  is  comprised  of  the  first  several  Fresnel 
zones  [52]  and  is  typically  small.  Image  theory  can  be  used  to  determine  the  contribution  due  to  the 
specu  ar  point.  Surrounding  the  specular  point  is  the  glistening  surface  that  is  comprised  of  quasi- 
specular  points  formed  by  favorably  tilted  planes  or  facets.  This  is  an  area  of  strong  scattering  present  for 
small  out-of^plane  angles  near  zero.  The  scattering  mechanism  at  work  in  the  glistening  surface  is  plane 
boundary  reflections  from  tilted  planes  that  are  physically  large  compared  to  a  wavelength.  (Most  of  us 
ave  expenenced  an  optical  version  of  a  glistening  surface  when  viewing  the  reflection  of  the  sun  or 
moon  in  a  body  of  slightly  rough  water.)  The  third  area  of  bistatic  scattering  is  the  diffuse  scattering 
which  IS  the  weakest  form  of  scattering.  Here  the  out-of-plane  angle  is  large  so  that,  neglecting  discrete 
sea  erers,  t  ere  are  usually  no  planes  favorably  tilted  for  generating  quasi-specular  returns.  Instead,  the 
scattenng  mechanism  is  Bragg  resonant  scattering  which  is  characterized  by  small  perturbation  theory 
^d  the  scattenng  strength  is  proportional  to  the  small-scale  surface  roughness  spectrum.  Here  the 
deviations  of  the  surface  height  about  its  mean  are  assumed  to  be  smaller  than  a  wavelength.  More  detail 
on  bistatic  scattenng  and  bistatic  clutter  modeling  can  be  found  elsewhere  [10],[53]-[55].  Note  that  the 
ocahzation  techniques  that  are  the  subject  of  this  report  do  not  rely  on  a  specific  scattering  model  or 
mechanism;  the  only  requirement  is  that  out-of-plane  bistatic  scattering  is  present. 
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The  out-of-plane  scattering  is  the  weakest  of  the  three  types  of  bistatic  scattering  described  above, 
but  it  is  this  weak  scattering  that  will  be  used  for  localization.  A  principal  reason  why  the  out-of-plane 
scattering  is  useful  is  that  the  scattered  signals  are  easily  resolved,  in  delay,  Doppler,  and  azimuth,  from 
the  strong  direct-path  signal.  It  is  generally  not  possible  to  resolve  the  glistening  surface  or  the  specular 
multipath  point  from  the  direct  path  because  the  differential  time  of  arrival  is  too  small  to  be  observed 
with  narrowband  receivers.  For  example,  using  the  altitudes  for  Mountaintop  in  Table  2.1  and  an  emitter 
slant  range  of  200  km,  the  path  length  of  the  in-plane  specular  multipath  route  is  200.096  km.  The 
difference  of  96  m  corresponds  to  a  time  difference  of  arrival  of  only  0.32  ps.  This  would  require 
receiver  bandwidths  exceeding  3  MHz  to  resolve  using  standard  correlation  techniques.  A  similar 
argument  holds  for  angle  measurements.  The  specular  multipath  and  glistening  surface  are  at  the  same 
azimuth  angle  as  the  direct  path  and  the  separation  in  elevation  is  a  small  fraction  of  the  elevation 
beamwidth  that  cannot  be  resolved  using  conventional  beamforming  techniques. 


2.3.1  Direct  and  Scattered  Signal  Levels 


The  signal  transmitted  by  the  emitter  travels  to  the  receiver  via  the  direct  propagation  path  and 
reflections  from  the  ground.  The  received  power  due  to  the  direct  path  is 


P,G,G^^ 

{AnRf 


(2.1) 


where  P,  is  the  emitter  transmit  power,  Gt  is  the  transmit  antenna  gain,  Gr  is  the  receive  antenna  gain,  R  is 
the  range  from  the  emitter  to  the  receiver,  and  A  is  the  receiver  wavelength.  The  aperture  of  the  receive 
antenna  Ar  is  related  to  the  gain  by  Ar  =  Gr^^Mn,  The  effective  radiated  power  (ERP)  is  the  product  of  the 
transmitter  power  Pt  and  the  transmitter  gain  Gf. 


The  average  power  reflected  from  a  differential  patch  of  ground  dPr  is  proportional  to  the 
scattering  cross  section  per  unit  area,  Ob,  times  the  area  of  the  patch,  ds.  The  total  average  power  received 
is  obtained  by  integrating  the  differential  power  over  the  area  illuminated  by  the  transmitter  and  receiver 
antenna  patterns  in 


P  = 


(2.2) 


where  R,  and  Rr  are,  respectively,  the  ranges  from  the  transmitter  and  the  receiver  to  the  differential  patch 
of  ground  ds.  This  expression  is  the  area-extensive  bistatic  radar  equation  [47].  The  received  power  from 
a  discrete  scatterer,  such  as  a  water  tower,  building,  or  terrain  feature,  can  be  obtained  using  the  standard 
bistatic  radar  equation  given  by 


P.  = 


P,G,G^X^<y 
(AnfRfR^  ' 


(2.3) 
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where  <T  is  the  radar  cross  section  (RCS)  of  the  discrete  scatterer.  The  discrete  scatterers  are  commonly 
modeled  as  uniformly  distributed  scatterers  with  densities  of  one  40-dBsm  scatterer  per  square  mile  [56] 
for  higher  frequencies  (X-band).  Results  described  in  Section  6.3  suggest  discrete  clutter  distributions 
closer  to  0.01  scatterers  per  square  mile  for  the  Mountaintop  system  and  parameters  at  WSMR. 

We  now  consider  the  relative  amplitudes  of  the  signals  described  by  (2.1)-(2.3)  using  the 
Mountaintop  system  parameters  in  Table  2.1.  An  emitter  with  an  ERP  of  2  kW  located  at  a  range  of  200 
km  will  result  in  a  received  direct-path  signal  of  —69  dBW  (decibels-Watt)  ignoring  the  specular  in-plane 
multipath.  The  receiver  thermal  noise  power  is  equal  to  kToBF^,  where  k  is  Boltzmann’s  constant,  To  is 
the  temperature,  B  is  the  receiver  bandwidth,  and  F„  is  the  receiver  noise  figure.  (At  room  temperature, 
^^0  =  4x10^*  W/Hz.)  The  resulting  noise  power  is  —143  dBW  and  hence  the  signal-to-noise  ratio  (SNR) 
for  the  direct  path  is  74  dB  which  is  a  very  strong  signal.  Next  consider  an  application  of  (2.2)  to 
compute  the  power  received  from  a  100-km^  patch  of  ground.  The  average  reflectivity  observed  at 
WSMR  of  -45  dB  [6], [9]  will  be  assumed.  The  RCS  cr  is  equal  to  the  product  of  the  reflectivity,  Cg,  and 
the  area  (80  dBm^)  and  is  equal  to  35  dBm^.  If  we  assume  the  ranges  R,  =  Rr  =150  km,  then  (2.3)  yields  a 
received  power  of  —144.5  dBW  for  an  SNR  of  —1.5  dB  which  is  a  weak  signal.  Next  we  consider 
scattering  from  a  clutter  discrete  at  the  same  range  as  the  clutter  patch.  Discrete  scatterers  with  cross 
sections  as  high  as  60  dBm^have  been  observed  [9],  but  we  shall  assume  a  more  conservative  level  of  40 
dBm^.  Equation  (2.3)  yields  a  received  power  level  of  -139.5  dBW  for  an  SNR  of  3.5  dB  which  is  also  a 
weak  signal.  Signal  integration  techniques  will  be  used  to  boost  these  signal-to-noise  ratios  to  values  that 
will  be  more  useful.  For  example,  if  1000  samples  can  be  coherently  integrated,  then  these  ratios  will  be 
boosted  by  30  dB  (assuming  an  integration  gain  equal  to  the  number  of  samples,  N). 

2.4  LOCALIZATION  OF  A  MOVING  EMITTER 

The  remainder  of  this  section  contains  a  description  of  the  bistatic  delay  and  Doppler  contours  and 
the  delay-Doppler  locus  or  ridge  that  is  observed  in  a  multipath  receive  beam.  The  bistatic  Hough 
transform  is  introduced  as  a  technique  for  localizing  an  emitter  under  an  idealized  scattering  scenario. 
This  technique  is  useful  because  it  provides  insight  by  permitting  a  graphical  interpretation  to  the 
localization  process.  More  rigorous  statistical  formulations  will  be  developed  in  Chapters  3  and  4. 

2.4.1  Delay-Doppler  Ridge 

Differential  delay  and  Doppler  contours  for  a  moving  emitter  and  a  stationary  receiver  are  shown 
in  Figure  2.3.  The  contours  of  constant  Doppler  are  hyperbolas  formed  by  the  intersections  of  cones  with 
the  ground.  We  shall  assume  that  the  ground  is  flat.  The  cones  are  described  by  f  =  V jX  cosy/ ,  where 
V  is  the  emitter  velocity,  A  is  the  emitter  wavelength,  and  y/  is  the  angle  measured  from  the  emitter 
direction  vector.  The  contoms  of  constant  delay  are  intersections  of  ellipsoids  with  the  ground.  If  both 
the  receiver  and  the  emitter  are  located  on  the  ground,  then  the  contours  become  ellipses  with  foci  at  their 
locations. 
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Kilometer 


Figure  2.3.  Bistatic  delay  and  Doppler  contours  for  an  emitter  at  a  range  of  150  km  traveling  north  at  a  velocity 
corresponding  to  a  maximum  Doppler  shift  of  200  Hz  directly  in  front  of  the  emitter  and  -200  Hz  behind  the 
emitter.  The  Doppler  shift  for  the  direct-path  signal  is  zero  for  this  geometry  and  the  heading  is  90  degrees. 


The  receiver  depicted  in  Figure  2.3  has  formed  a  beam  in  the  direction  <!>  which  is  30  degrees.  This 
beam  is  assumed  to  be  very  narrow  and  is  shown  in  the  figme  as  a  line.  The  emitter  signal  scattered  by  the 
ground  at  each  point  along  the  line  will  be  delayed  and  Doppler  shifted  by  amounts  indicated  by  the 
differential  delay  and  Doppler  contours.  The  energy  in  the  scattered  signal  will  be  confined  to  a  locus  or 
ridge  in  delay-Dopplar  space  as  shown  in  Figure  2.4.  The  remainder  of  the  scattered  delay-Doppler 
intensity  map  will  contain  only  receiver  thermal  noise.  The  diffo-ential  Doppler  ranges  from  zero  Hz 
corresponding  to  the  direct  path,  to  a  maximum  of  200  Hz  at  a  bistatic  dday  near  100  km  The  ridge  is 
shown  as  a  line  since  the  antenna  beamwidth  was  assumed  to  be  very  narrow.  These  delay-Doppler  ridges 
have  been  observed  in  practice  [9]  whare  the  ridge  was  spread  somewhat  due  to  the  nonzero  antenna 
beamwidths.  In  the  next  section  we  use  the  bistatic  Hough  transform  to  localize  an  emitter  by  processing 
several  scatterers  on  the  delay-Doppler  ridge. 
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Figure  2.4.  Delay-Doppler  ridge  for  bistatic  scattering  from  a  moving  emitter  and  the  geometry  of  Figure  2.3. 


2.4.2  Bistatic  Hough  Transform 

In  this  section  we  derive  the  bistatic  Hough  transform.  This  is  a  numerical  transform  technique  that 
is  closely  related  to  template  matching  [57]  and  the  Radon  transform  [58].  It  will  be  used  to  localize  an 
emitter  in  an  idealized  scattering  scenario  and  to  provide  insight  into  the  loceilization  using  out-of-plane 
multipath  process.  It  will  also  be  useful  as  a  tool  for  exanoining  ambiguities  that  may  occur  when  using 
out-of-plane  multipath  for  emitter  localization. 

The  Hough  transform  (HT)  was  introduced  as  a  means  to  detect  straight  line  segments  in  binary 
digital  images  [59].  It  maps  a  spatially  extended  shape  such  as  a  line  in  image  space  to  a  point  in 
parameter  space.  The  first  application  of  the  HT  for  use  in  radar  was  considered  in  [60]  where  it  was 
used  to  detect  weak  linear  target  tracks.  In  this  section  the  HT  will  be  used  to  detect  the  bistatic  delay- 
Doppler  ridge. 

The  Hough  transform  for  line  detection  can  be  written  in  integral  form  as  [61] 
H(p,0)=jjl(x,y)8(p-xcosG -ysm9)dxdy  ,  (2.4) 

where  I(x,  y)  represents  the  image  space  and  H(p,0}  represents  the  parameter  space.  The  impulse  function 
is  nonzero  only  when  p  =  x  cos  6  +  y  sin  0 .  This  expression  in  known  as  the  normal  parameterization  of 
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a  line.  This  parameterization  is  preferred  over  the  more  common  parameterization,  y  =  mx  +  b,  because  it 
avoids  the  numerical  problems  that  can  arise  when  dealing  with  vertical  lines  which  have  infinite  slopes. 
For  a  binary  image  I(x,  y)  and  the  f  pair  of  parameters  (ptA),  the  integral  evaluates  to  the  number  of 
image  points  that  coincide  with  the  line  p^  =  xcosdi  +  ysin0,. .  For  example,  suppose  I(x,  yj  is  a  100  by 
100  binary  image  consisting  of  a  vertical  line  (0  =  7t/2)  passing  through  the  origin  (p=0).  The  Hough 
transform  will  contain  a  peak  in  parameter  space  of  magnitude  100  located  at  —  (0,^/2)  and 

some  very  small  residual  terms — mostly  equal  to  one  or  zero  at  other  values  of  (/j,0)  . 


For  the  emitter  localization  problem,  the  image  space  is  the  delay-Doppler  map  and  the  parameter 
space  contains  the  emitter  parameters  we  want  to  estimate.  A  bistatic  Hough  transform  designed  to  detect 
the  delay-Doppler  ridge  contained  in  a  bistatic  delay-Doppler  map  can  be  written  as 

/  A 


R- 


_ D _ 

^iny(g,y,/)  ^  sin<^ 
^sma{H,V,f)  sma{H,V,f)  ^ 


\iDdf  , 


) 


(2.5) 


where:  H{R,  H,  V)  is  the  HT  in  parameters  R,  H,  V, 

M(D,f)  is  the  bistatic  delay-Doppler  intensity  map, 

D  is  the  differential  range  delay, 

/is  the  differential  Doppler  shift  (measured  relative  the  direct  path), 

R  is  the  range  from  the  receiver  to  the  emitter, 

H  is  the  emitter  heading, 

V  is  the  emitter  velocity, 

(p  is  the  radar  look  azimuth, 

7  is  the  angle  to  the  scatterer  which  is  a  function  of  D,  H,  and  V,  and 

cc  is  the  remaining  angle  in  the  bistatic  triangle  which  is  shown  in  Figure  2.5. 


Scatterer 


Figure  2.5.  Bistatic  triangle  formed  by  the  receiver,  emitter,  and  scatterer. 
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The  argument  of  the  impulse  function  is  taken  from  (1.2).  This  equation  was  derived  using  the  law 
of  sines  to  compute  the  emitter  range  using  a  single  bistatic  scatterer  when  the  emitter  velocity  and 
heading  is  known.  Here  the  emitter  velocity  and  heading  are  parameters  that  will  be  treated  as 
independent  variables  for  the  integration  in  (2.5).  The  integrand  is  non-zero  only  when  the  argument  of 
the  impulse  function  is  equal  to  zero.  To  evaluate  the  integral,  we  solve  for  as  a  function  of  the 
measured  quantities  of  delay  Z),  differential  Doppler  /,  and  radar  azimuth  0;  and  the  adjustable 
parameters  of  emitter  velocity  V  and  heading  H.  The  angle  yis  related  to  the  emitter  velocity,  heading, 
and  differential  Doppler  and  can  be  determined  as  follows. 

The  direct-path  absolute  Doppler  shift  is 

fd=jCOs(m-H)  ,  (2.6) 

emd  the  scatterer  absolute  Doppler  is 

V  V 

fs=jCos(y/)  =  jCos(lSO-H-y)  .  (2.7) 

The  differential  Doppler  that  is  measured  by  the  receiver  is 

/  =  /^-/d  =— [cos(180-//-7)-cos(180-//)]  ,  (2.8) 

but  cos(180-/f)  =  -cos(if)  and  cos(180- = -cos(//  +  7)  ,  so 

/  =  y[-cos(/t+7)  +  cos(H)]  .  (2.9) 


Equation  (2.9)  can  be  used  to  determine  /since 

^cos(H  +  y)  =  ^cos(H)-f  , 


(2.10) 


and  finally 


7  =  -H+cos  *  cos(H)  - 


(2.11) 


a  =  180-^- 7  . 


(2.12) 


Using  the  above  equations  (2.5)— (2.12),  we  can  describe  the  Hough  transform  algorithm  with  pseudo¬ 
code  in  Figure  2.6. 
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For  each  delay-Doppler  cell  M(D,f)  containing  scattered  emitter  energy 
For  each  trial  emitter  velocity  parameter  V 
For  each  trial  emitter  heading  H 

calculate  R  (from  argument  of  impulse  function  in  (2.5) 
Accumulate  counts  in  parameter  space  {R,  H,  V) 

end 

end 

end 

Search  parameter  space  (R,H,V)  for  peak. 

Figure  2.6.  Pseudo -code  for  bistatic  Hough  transform  localization  algorithm. 


Localization  using  the  Hough  transform  will  now  be  demonstrated  by  an  example.  The  emitter  is 
located  at  a  range  of  150  km  traveling  with  a  heading  of  90  degrees  and  a  velocity  of  140  m/s  which 
corresponds  to  a  maximum  absolute  Doppler  shift  of  200  Hz.  We  first  assume  that  the  ermtter  velocity  is 
known.  This  restricts  the  parameter  space  to  only  two  dimensions  (R,  H)  thereby  reducing  the  parameter 
search  to  only  two  dimensions. 

Consider  first  a  single  scatterer  located  on  the  delay-Doppler  ridge  as  shown  in  Figure  2.7(a).  The 
Hough  transform  of  this  single  point  in  delay-Doppler  image  space  is  a  curved  line  in  parameter  space, 
which  is  shown  in  Figure  2.7(b).  This  result  may  be  interpreted  as  follows;  An  emitter  traveling  at  a 
known  velocity  cannot  be  localized  using  delay-Doppler  information  from  a  single  scatterer;  however, 
the  information  from  the  single  scatterer  confines  the  emitter  parameters  of  range  and  heading  to  the 
curve  in  Figure  2.7(b).  Transformation  of  additional  scatterers  in  delay-Doppler  space  will  result  in 
additional  curves  in  parameter  space  and  the  intersection  of  these  curves  yield  the  desired  emitter 
parameters.  This  is  illustrated  in  Figure  2.8  for  three  scatterers.  The  curves  intersect  at  the  correct  emitter 
range  of  150  km  and  heading  of  90  degrees.  This  is  somewhat  similar  to  the  intersecting  hyperboloids 
discussed  in  Section  1.1.1. 

When  many  scatterers  are  present,  a  sharp  peak  occurs  in  the  parameter  space.  This  is  shown  in 
Figure  2.9  for  about  30  scatterers.  The  emitter  parameters  are  estimated  by  searehing  the  parameter  space 
for  the  peak.  For  this  example,  the  impulsive  peak  clearly  identifies  the  correct  emitter  range  of  150  km 
and  bearing  of  90  degrees. 
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Differential  Doppler  (Hz) 


(a)  (b) 


Figure  2. 7.  Hough  transform  of  (a)  point  in  delay-Doppler  space  is  (b)a  curve  in  the  trial  range-heading  space. 


(a)  (b) 


Figure  2.8.  (a)  Three  scatterers  located  on  the  delay-Doppler  locus  and  (b)  the  corresponding  Hough  transform. 


Trial  Range  (km) 

(c)  (d) 


Figure  2.9.  Localization  using  Hough  transform,  (a)  Thirty  scatterers  on  the  delay-Doppler  ridge  are  transformed 
to  emitter  parameter  space  (b,  c,  d).  Parts  (c)  and  (d)  are  views  of  the  parameter  space  after  the  curves  shown  in  (b) 
have  been  accumulated. 
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Figure  2.9(b)  can  be  studied  from  a  geometric  dilution  of  precision  (GDOP)  point  of  view.  GDOP 
is  defined  as  the  ratio  of  the  RMS  position  errors,  e.g.,  emitter  range,  and  the  measurement  errors,  e.g., 
differential  delay  (in  km)  [62].  A  large  GDOP  means  that  small  measurement  errors  will  produce  large 
parameter  errors.  For  example,  consider  the  HT  of  the  two  scatterers  at  the  longest  delay.  This  consists  of 
the  two  most  vertical  lines  in  parameter  space  which  are  seen  to  be  nearly  parallel  in  Figure  2.9(b).  While 
these  lines  intersect  at  the  correct  point  in  parameter  space,  it  is  easy  to  see  that  a  small  measurement 
error  in  image  space  resulting  in  the  movement  of  one  of  these  lines  will  produce  a  large  error  in  the 
parameter  space.  Similarly,  the  HT  of  the  two  shortest  delay  scatterers  are  also  parallel  curves.  However, 
if  we  consider  the  lines  due  to  a  short-delay  scatter  and  a  long  delay  scatterer,  we  can  see  that  these  are 
nearly  perpendicular  and  thus  small  measurement  errors  will  produce  small  parameter  errors.  The  effects 
of  measurement  errors  on  parameter  estimation  errors  will  be  examined  in  Chapter  5. 

2.4.3  Uniqueness  of  Parameter  Space  Mapping 

In  the  previous  section  an  emitter  was  localized  in  the  two  dimensions  of  range  and  heading  where 
the  velocity  was  assumed  to  be  known.  The  HT  was  used  to  perform  the  localization  which  generated  a 
single,  unambiguous  result  using  only  one  beam  of  received  scattering  data.  We  now  consider  the  case 
where  all  three  parameters  are  unknown  and  show  that  multiple  beams  may  be  required  to  generate 
unique,  unambiguous  localization  results. 

The  HT  was  determined  for  three  trial  emitter  velocities  corresponding  to  maximum  Doppler  shifts 
of  160,  200,  and  240  Hz.  The  value  of  200  Hz  corresponds  to  the  true  emitter  velocity.  We  expect  a  peak 
at  the  true  velocity  corresponding  to  200  Hz  and  small  residual  terms  for  the  two  other  trial  velocities. 
The  HT  for  the  three  cases  is  shown  in  Figure  2.10.  We  see  that  for  the  true  velocity,  the  lines  converge 
at  a  single  point  that  will  generate  a  sharp  peak  in  parameter  space.  For  the  other  two  velocities,  the  lines 
do  not  converge  as  sharply,  suggesting  blurred  or  smeared  peaks  in  parameter  space.  For  this  noise-free 
case,  we  can  cMm  that  the  HT  has  successfully  localized  this  emitter  in  three  dimensions  of  range, 
heading,  and  velocity;  however,  if  noise  was  present,  the  true  peak  could  be  blurred  as  well.  Hence  it  may 
be  difficult  to  select  the  correct  emitter  parameter  set. 

To  see  why  peaks  in  parameter  space  occurred  for  three  different  parameter  sets,  we  generate 
delay-Doppler  ridges  corresponding  to  each  one.  The  false  parameter  sets  are  R=1 16  km,  H=76  deg.,  and 
V=80%  of  the  true  velocity,  and  are  R=177  km,  H=99.5  deg.,  and  V=120%  of  the  true  velocity.  The 
delay-Doppler  ridge  for  the  three  parameter  sets  are  shown  in  Figure  2.11.  We  can  see  that  the  three 
ridges  are  nearly  identical. 
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Doppler  (Hz) 


Trial  Heading  (deg) 

Figure  2.10.  Hough  transforms  of  a  measured  delay-Doppler  ridge  assuming  three  different  emitter  velocities. 


Delay  (km) 


Figure  2.11.  Delay-Doppler  ridges  for  the  three  parameter  sets  suggested  by  Figure  2.10.  The  delay-Doppler  ridges 
are  nearly  identical. 


This  result  suggests  that  the  localization  of  this  emitter  may  be  difficult  because  we  have  shown 
that  three  different  parameter  sets  possess  nearly  identical  delay-Doppler  ridges.  This  is  an  important 
result  because  the  localization  technique  relies  on  the  uniqueness  of  the  delay-Doppler  ridge  to  estimate 
the  emitter  parameters.  However,  up  to  this  point  we  have  considered  only  the  contributions  from 
scattering  in  a  single  receive  beam  at  30  degrees.  When  additional  beams  are  employed,  the  mismatch 
between  the  true  parameters  and  the  false  peaks  will  be  increased.  For  example,  if  we  form  an  additional 
radar  beam,  say  at  -30  degrees,  we  see  in  Figure  2.12  that  the  delay-Doppler  ridges  for  these  three 
parameter  sets  are  no  longer  similar.  Hence  only  the  true  peak  will  be  enhanced  by  the  additional  beam 
and  the  localization  result  will  be  unique  as  expected. 


Figure  2.12.  Delay-Doppler  ridges  for  the  three  parameter  sets  and  a  —30  degree  receive  azimuth. 


2.5  CHAPTER  SUMMARY 

This  section  contains  the  background  and  introductory  information  required  to  develop  the 
estimation  techniques  described  in  the  following  chapters.  A  discussion  of  signal  and  noise  models  and 
bistatic  scattering  is  provided.  The  delay-Doppler  ridge  is  described  and  the  bistatic  Hough  transform  is 
introduced  and  used  to  localize  an  emitter  in  an  idealized  case.  This  technique  is  graphical  in  nature  and 
promotes  a  better  understanding  of  the  localization  process.  Statistically  based  solutions  are  provided  in 
subsequent  chapters.  The  HT  is  also  shown  to  be  useful  as  a  tool  for  examining  potential  localization 
ambiguities.  We  used  it  to  show  that  multiple  receive  beams  may  be  required  to  perform  emitter 
localization. 
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3.  LOCALIZATION  IN  HOMOGENEOUS  CLUTTER 


3.1  INTRODUCTION 

In  this  chapter  we  design  an  estimator  based  on  an  assumption  that  the  bistatic  clutter  is 
homogeneous.  The  homogeneous  clutter  model  will  be  described  along  with  several  assumptions  that 
lead  to  a  simple  statistical  formulation.  The  likelihood  function  and  corresponding  maximum  likelihood 
estimator  (MLE)  will  be  determined.  The  structure  of  the  estimator  is  shown  to  be  a  time-varying  FIR 
filter  with  Doppler  modulators  chosen  to  match  the  Doppler  shifts  in  the  multipath  signals.  The  optimum 
solution  technique  will  be  described  as  well  as  a  sub-optimum  technique  that  employs  weight  thinning. 


3.2  BASIC  SCATTERING  MODEL 


The  fundamentals  of  bistatic  scattering  were  introduced  in  Section  2.3.  For  the  homogeneous 
scattering  model  assumption,  the  average  received  power  for  a  delay-angle  cell  was  written  as 


area 


P,G,G,k^a„ 


■ds  , 


(2.2) 


where  the  reflectivity  Oo  does  not  vary  much  from  cell  to  cell.  This  describes  area  scattering  because  the 
received  power  depends  on  the  product  of  the  reflectivity  and  the  area  of  the  delay-azimuth  cell.  This 
differs  from  point  scattering  which  originates  from  a  discrete  scatter  with  RCS  ct.  For  the  homogeneous 
scattering  assumption,  we  assume  that  the  area  scattering  is  stronger  than  the  discrete  scattering.  In  the 
following  chapter,  we  treat  the  inhomogeneous  case  where  we  assume  that  point  scattering  is  stronger 
than  area  scattering.  These  two  assumptions  are  proposed  because  they  permit  fairly  straightforward 
estimator  designs.  In  practice,  bistatic  clutter  will  exhibit  properties  of  both  area  and  point  scattering.  The 
preferred  estimator  may  therefore  depend  on  the  scattering  characteristics  of  the  specific  terrain  for  any 
given  scenario. 

The  bistatic  reflectivity  Oq  in  (2.2)  is  generally  a  function  of  the  incidence  angles  and  the  out-of¬ 
plane  angle.  In  the  out-of-plane  region  the  simplest  bistatic  clutter  models  are  based  on  the  constant- 
gamma  model  for  monostatic  clutter  [63] 

(T^  =7sin0;  ,  (3.1) 

where  <7q  is  the  monostatic  cross  section  per  unit  area,  0,  is  the  incidence  angle,  and  yis  an  additional 
scattering  parameter  that  is  assumed  to  be  independent  of  0,.  This  may  be  generalized  for  bistatic 
scattering  [63]  as 

(Tq  =  7(sin0^  sin0,y'^^  ,  (3-2) 
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where  Or  and  0,  are  the  receiver  incidence  angle  and  transmitter  incidence  angle,  respectively.  Values  of 
the  reflectivity  parameter  7  were  found  to  be  approximately  -20  dB  for  the  Mountaintop  system  at 
WSMR  [9]. 

3.3  STATISTICAL  FORMULATION 

The  statistical  formulation  begins  with  a  statement  of  the  clutter  model  and  ends  with  the 
likelihood  function.  In  order  to  write  the  likelihood  function,  several  further  simplifications  and 
approximations  are  required. 

We  shall  assume  that  the  antenna  beamwidth  is  very  narrow  and  the  delay  cells  are  closely  spaced. 
This  helps  to  justify  the  further  assumption  that  the  scattering  from  each  delay-azimuth  cell  originates 
from  the  center  of  the  cell.  This  permits  us  to  ignore  Doppler  dispersion  (broadening)  and  is  illustrated  in 
the  delay-Doppler  map  shown  in  Figure  3.1  where  the  scattering  centers  in  the  multipath  beam  are 
denoted  by  dots.  The  significance  of  this  assumption  is  that  it  permits  us  to  treat  the  differential  delay  and 
azimuth  for  each  scattering  cell  as  independent  variables.  The  delay  is  simply  the  center  of  each  delay  bin 
while  the  azimuth  is  the  look  direction  of  the  beam.  (For  the  inhomogeneous  scattering  model  described 
in  the  next  section,  this  assumption  will  not  be  valid.  There  we  will  be  required  to  estimate  the  delay  and 
azimuth  for  each  discrete  scatterer  separately.) 

The  very  high  signal-to-noise  ratio  of  the  emitter  waveform  via  the  direct  propagation  path  permits 
additional  simplifications.  In  Section  2.3.1  we  saw  that  the  SNR  of  the  direct  path  was  approximately  74 
dB  for  the  typical  Mountaintop  system  parameters.  The  direct-path  signal  is  measured  by  the  reference 
beam  which  is  pointed  directly  at  the  emitter.  The  SNR  of  the  emitter  signal  that  was  scattered  by  a 
100-km^  patch  of  ground  was  approximately  -1.5  dB  per  pulse.  The  very  high  SNR  for  the  emitter  signal 
via  the  direct  path  facilitates  a  very  accurate  sampling  of  the  emitter  signal.  For  the  purposes  of  this 
section,  we  will  consider  this  signal  to  be  known  exactly.  Furthermore,  the  very  high  emitter  SNR 
facilitates  emitter  azimuth  estimates  that  are  very  accurate.  Again,  for  the  purposes  of  this  section,  we 
will  assume  that  azimuth  of  the  emitter  is  known  exactly. 

With  the  above  assumptions  and  simplifications,  we  are  ready  to  write  an  expression  for  the 
received  signal  in  the  multipath  beam.  We  write 

b  =  \,2,3,..,B  ,  (3.3) 

p=\ 

where  Zb(n)  is  the  n  sample  from  the  beam,  Ubp  and  co^p,  respectively,  represent  the  complex 
amplitude  and  Doppler  frequency  in  radians  per  second  of  the  independent  scattering  cell  for  the 
independent  beam.  The  differential  delay  Dp  here  is  an  integer  equal  to  p-\.  The  emitter  waveform  is  s(n) 
and  the  zero-mean  complex  observation  noise  is  denoted  by  gb(n}  with  covariance  matrix  Rgg.  The 
dimension  of  Rgg  Nhy  N  where  N  is  the  total  number  of  samples  in  the  data  vector  Zb.  In  general,  the 
noise  term  gb(n)  contains  anything  not  included  in  the  summation  term  representing  the  multipath  model. 
For  example,  scattered  signals  that  enter  via  the  antenna  sidelobes,  distant  scatter  from  beyond  the 


30 


summation  length,  additional  emitter  signals,  etc.,  are  all  included  in  the  noise  term.  These  noise  sources 
are  referred  to  as  model  noise  because  they  are  unwanted  signal  soinces  not  described  by  the  model.  For 
our  purposes  in  this  section,  we  shall  assume  that  the  noise  is  receiver  thermal  noise  which  is  white  and 
uncorrelated  from  beam-to-beam  and  from  sair^le-to-san5)le. 
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Figure  3.1.  Range  and  Doppler  contours  for  a  moving  emitter.  A  reference  beam  pointed  at  the  emitter  direct 
samples  the  emitter  waveform  and  a  multipath  beam  collects  the  attenuated,  delayed,  and  Doppler  shifted  replic 
of  the  emitter  waveform. 


We  can  write. (3. 3)  in  vector  form  as 

Zb=«»b+gb  -  (3.4) 

where  each  term  is  an  by  1  vector  and  mifn)  is 

b  =  l,2X.,B  .  (3.5) 

p=\ 


31 


That  is,  Zb  is  equal  to  mb  when  the  measurement  noise  gb  is  equal  to  zero.  The  probability  density 
function  for  the  measurements  Zb(n)  can  be  written  as  a  complex  Gaussian 


L 


b 


1 


(3.6) 


for  a  single  beam,  and  as 


4w=n 

b=l 


(3.7) 


or 


4(z)= 


b=l 


(3.8) 


for  multiple  beams  where  we  have  assumed  that  the  covariance  matrix  of  the  noise,  R,  is  the  same  for 
each  beam.  (The  gg  subscript  in  Rgg  has  been  dropped  for  convenience.)  For  the  white  thermal  noise 
assumption,  R  is  of  the  form  ,  and  we  can  write  the  PDF  of  the  measured  vector  for  a  given  beam 
(3.6)  as 


4(H^)= 


1 


yTt  a  j 


-(z-in(0))"  (z-m(0))/t7' 


(3.9) 


where  we  have  dropped  the  beam  subscripts.  This  may  also  be  written  as 


\  \  l(J^ 


71  G 


(3.10) 


where  |  |  denotes  the  vector  norm.  The  term  m(^  is  the  N  x  1  vector  whose  elements  for  a  given  beam 
are  given  in  (3.5)  by  setting  n  =  1,2, ..,N.  The  parameter  vector  6  has  been  included  to  remind  the  reader 
that  this  vector  is  a  function  of  the  emitter  parameters  we  want  to  estimate.  Conditional  joint  probability 
or  likelihood  functions  of  the  form  of  (3.9)  or  (3.10)  will  be  assumed  throughout  the  rest  of  this  chapter. 
Note  that  these  are  conditional  joint  probability  functions  when  the  parameters  0  are  fixed  and  the 
observation  vector  z  is  variable,  and  these  are  likelihood  functions  when  the  observation  vector  z  is  fixed 
2ind  parameters  0are  variable.  The  nature  of  the  dependence  of  (z|0)  on  6  is  in 


p=l 


(3.11) 


where  the  only  term  that  depends  on  0is  (Op(0),  the  Doppler  frequency  shift  for  each  scattering  cell.  That 
is  a  direct  consequence  of  the  homogeneous  clutter  assumption  which,  as  mentioned  above,  permits  us  to 
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treat  azimuth  and  time  delay  as  independent  variables.  We  still  need  to  estimate  the  complex  amplitude. 
Op,  for  each  scattering  cell.  The  scatterer  amplitudes  are  nuisance  parameters. 

3.4  MAXIMUM  LIKELIHOOD  ESTIMATOR 

In  this  section  we  determine  the  maximum  likelihood  estimator  (MLE)  of  the  parameter  vector  6, 
based  on  the  model  described  above.  The  MLE  produces  the  most  plausible  estimate  vector  6  given  a 
measurement  vector  z  by  maximizing  the  likelihood  function  (3.10).  Maximizing  this  function  is 
equivalent  to  minimizing  the  exponent 

|z-in(0)r  =  |gr  ,  (3-12) 

where  the  vector  iii(6)  is  the  system  model  from  (3.4)  and  (3.11)  and  z  is  the  measured  data  vector.  The 
difference  between  these  two  terms  is  the  observation  error,  g.  We  desire  to  determine  the  parameter 
vector  6,  as  well  as  the  scatterer  weights  Up,  that  yield  the  minimum  error.  We  will  be  required  to  search 
the  parameter  space  while  estimating  the  nuisance  parameters  Up  for  the  minimum  squared  error.  The 
parameter  set  that  minimizes  the  squared  error  will  be  chosen  as  the  localization  estimate. 

We  first  treat  the  case  where  the  Doppler  terms  cOp  are  zero  and  then  extend  the  results  to  handle 
nonzero  case.  When  the  Doppler  terms  are  zero,  (3.5)  becomes 

p 

m(n)=^apS(n  — Dp)  .  (3.13) 

p=i 

Hence  this  problem  reduces  to  the  standard  system  identification  problem  [64]  depicted  in  Figure  3.2. 
The  multipath  channel  is  modeled  by  a  discrete-time  linear,  time-invariant  FIR  filter  with  filter 
coefficients  Op.  The  input  is  the  emitter  waveform  s(n),  the  output  of  the  noisy  multipath  channel  is  z(n), 
and  the  system  identification  filter  output  is  m(n). 


Emitter 

signal 


s(n) 


Emitter 

Multipath 

Channel 

z(n) 

i 

.i 

FIR 

System 

Model 

m(n)  1 

9(n) 


Figure  3.2.  Finite  impulse  response  (FIR)  model  for  system  identification. 
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The  optimum  filter  coefficients  (in  a  least-squares  sense)  minimize  the  sum  of  the  square  error 
sequence  g(n), 


•^=2k(«)r 

n=P 

N  P 

=2 


(3.14) 


The  summation  over  n  starts  at  P  to  assure  that  all  the  P  filter  coefficients  have  nonzero  input  data 
values.  To  determine  the  filter  coefficients  that  minimize  J,  we  differentiate  (3.14)  with  respect  to  the  real 
and  imaginary  parts  of  each  filter  coefficient  Up  and  set  the  result  equal  to  zero.  The  optimum  filter 
coefficients  are  the  solutions  to 


N 


(3.15) 


which  is  the  temporal  version  of  the  principle  of  orthogonality.  The  expression  on  the  left  side  of  (3.15) 
represents  the  time  average  of  the  cross-correlation  between  the  p'*  tap  input  s{n  -  Dp)  and  the  minimum 
estimation  error  §„,„(«)  over  all  values  of  the  time  index  n  in  the  interval  (P,  N).  The  principle  of 
orthogonality  states  that  when  the  filter  coefficients  are  selected  to  minimize  J,  then  the  minimum  error 
series  {g„,„(n)}  is  orthogonal  to  input  time  series  {^(«  -  Dp) }  to  the  p'*  coefficient  of  the  filter  of  length  P. 
A  good  description  of  the  least-squares  problem  can  be  found  in  [65]. 


To  make  use  of  this  result,  we  rewrite  (3.3)  as 

p 

gin)=  zin)-  ^  flp  s(n  -  ) 

p=l 

and  substitute  for  g(n)  into  (3.15)  to  yield 

Y,^*in-DAz(n)-Y^a^s(n-D^)  =0,  A:  =  1,2,...,P  . 

n=P  p=l  J 

Expanding  and  rearranging  terms  yields 

n=l  n=P 


(3.16) 


(3.17) 


(3.18) 


where  the  inner  sum  on  the  left  is  the  time-average  autocorrelation  function  of  the  filter  input  for  lag  Dp- 
Dk  and  the  sum  on  the  right  is  the  time-average  cross-correlation  function  for  lag  A-  In  matrix  form,  we 
can  write  (3.18)  as 


R„a  =  r. 


(3.19) 
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(3.20) 


with  solution  for  the  optimum  weight  vector  given  by 

a  =  . 

which  is  the  well-known  result.  We  can  compute  and  using  the  (N~P +l)hy  P  data  matrix 


S(P) 

s(P-\)  - 

5(2)  5(1) 

s(P)  - 

5(2) 

siN-1) 

siN-2)  ••• 

: 

s(N) 

s(N-l)  ••• 

•••  s(N-P  +  l) 

as  =S"S ,  =S"z ,  and  a  =  (s"s)  ^S"z .  We  select  (N-P+l)>  P  so  that  S"S  is  invertable. 


The  above  development  generated  the  weight  vector  a  that  minimized  the  exponent  of  the 
likelihood  function  derived  for  Gaussian  noise.  Hence  these  are  maximum  likelihood  weights.  If  the 
distribution  of  the  observation  noise  was  non-Gaussian,  then  these  weights  are  optimum  in  the 
least-square  sense. 


How  can  the  above  results  be  extended  to  handle  the  time-varying  filter  suggested  by  (3.3)?  To 
solve  the  general  time-varying  case,  we  replace 


s{n-D  )hy  s{n- 


(3.22) 


in  (3.15)-(3.21).  Thus  we  have  added  the  differential  Doppler  term  which  contains  the  information 
necessary  to  estimate  the  parameter  vector  6  since  O  is  a  function  of  ft  The  actual  estimation  process 
will  require  a  search  over  all  possible  values  of  ft  in  parameter  space. 


When  we  perform  the  substitution  of  (3.22),  Equation  (3.17)  becomes 


N 

I 

n=P 


z{n)-Y.as(n-D)e 


=  ,  (3.23) 


p=\ 


and  (3.18)  becomes 


S 


jC0p(n-Dp) 


=^s*{n-D,)e  ^*^z(n)  fe=l,2,...,P  . 


Thus  the  individual  terms  in  the  autocorrelation  function  become 


(3.24) 
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(3.25) 


Rss^Kp)=I.s*{n-D^)e 

n=P 


,k=ia,...,p, 
P  =  12,...,P 


and 


^'‘^z(n),  k=l,2,...,P 


(3.26) 


and  finally,  the  optimum  weights 


a  =  R. 


(3.27) 


The  structure  of  this  estimator  is  the  time-vaiying,  FIR  tapped  delay  line  filter  that  is  illustrated  in 
Figure  3.3.  The  vector  a  represents  the  complex  filter  coefficients  or  taps  of  the  filter.  The  differential 
Doppler  terms  (Op  are  shown  modulating  the  delayed  replicas  of  the  emitter  waveform.  A  separate 
Doppler  modulator  cOp  is  employed  for  each  tap  Up. 


Figure  3.3.  Emitter  localization  filter  structure  for  a  single  multipath  beam.  The  emitter  parameter  set  that 
minimizes  the  output  power  is  the  ML  estimate. 


Note  that  we  are  not  trying  to  estimate  the  Doppler  terms  (Dp  directly.  Instead,  we  are  computing  a 
deterministic  set  of  modulators  fi)  for  each  trial  parameter  set  ft  We  then  compute  the  optimum  weights  a 
for  the  deterministic  modulators  and  evaluate  the  error  given  by  (3.14).  Since  the  Doppler  modulators 
will  contain  some  mismatch  for  all  cases  except  when  the  trial  emitter  parameters  equal  the  true  emitter 
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parameters,  the  error  can  be  large  even  though  the  weights  are  “optimum.”  When  the  trial  emitter 
parameters  equal  the  true  parameters,  the  error  is  minimized.  The  minimum  error  maximizes  likelihood 
function  in  (3.10)  and  is  hence  the  maximum  likelihood  estimator  for  the  emitter  parameters 

e={R,H,V). 

In  Section  2.4.3  we  saw  that  multiple  beams  may  be  required  to  resolve  potential  ambiguities.  The 
MLE  described  above  can  be  easily  extended  for  multiple  beams  as  shown  in  Figure  3.4.  This 
architecture  assumes  that  the  individual  multipath  beams  are  uncorrelated  (i.e.,  the  beams  are 


Emitter 


orthogonal).  Thus  we  will  be  required  to  solve  B  smaller  problems  of  size  P  instead  of  the  larger  problem 
of  size  BxP,  where  B  is  the  number  of  beams  and  P  is  the  number  of  taps. 

The  single  beam  structure  of  Figure  3.3  has  been  used  in  the  past  [9]  for  jammer  multipath 
mitigation  where  the  Doppler  modulators  were  used  to  improve  mitigation  performance.  There  the  error 
vector  was  the  target  beam  with  the  emitter  multipath  energy  removed.  It  is  re-introduced  here  as  the 
MLE  for  emitter  localization.  The  fundamental  difference  is  that  in  [9]  the  Doppler  modulators  were 
estimated  independently  and  used  to  improve  jammer  multipath  mitigation  performance.  For  localization, 
the  modulators  are  computed  separately  for  each  parameter  set  and  the  set  with  the  lowest  error 
(corresponding  to  the  best  mitigation)  is  selected  as  the  MLE.  Thus  we  see  that  jammer  multipath 
mitigation  and  emitter  localization  are  very  closely  related. 

The  only  thing  left  to  specify  in  this  estimator  is  how  to  compute  the  Doppler  modulators.  These 
are  given  by 


(o(R,H,V)= 


2kV 


( 


cos  H  —  cos 


H  -1-  sin' 


(Dl+2RDp)sm<j) 

Dl+2[R^+RDjl-cos(t>) 


(3.28) 


which  is  a  function  of  the  emitter  parameters  (R,  H,  V),  the  differential  delay  for  each  filter  tap  Dp,  the 
angle  of  the  multipath  beam  (j),  and  the  emitter  wavelength  A.  This  was  derived  based  on  a 
two-dimensional  approximation  that  ignores  the  altitude  of  the  emitter:  the  details  are  contained  in 
Appendix  A.  An  expression  based  on  a  three-dimensional  approximation  is  derived  in  Appendix  B  and 
the  exact  three-dimensional  solution  is  contained  in  Appendix  C. 


3.5  SOLUTION  TECHNIQUE 

The  fundamental  method  of  solution  is  the  dense  multidimensional  search  over  the  emitter 
parameter  space.  Hence  we  are  solving  this  inverse  problem  by  solving  an  “infinite”  number  of  forward 
solutions  and  selecting  the  one  that  maximizes  the  likelihood  function.  Such  dense  searches  are  very 
costly  especially  for  three  or  four  dimensions.  In  practice,  we  shall  first  perform  a  coarse  three- 
dimensional  search  over  emitter  range,  heading,  and  velocity.  The  peak  of  the  coarse  search  is  used  as  the 
starting  point  for  an  iterative  procedure,  such  as  the  Newton-Raphson  method.  The  iterative  procedure  is 
carried  out  in  all  four  dimensions  of  range,  heading,  velocity,  and  altitude.  Generally,  if  the  starting  point 
is  close  to  the  global  maximum,  the  iteration  will  converge  to  it  [67].  Thus,  the  initial  coarse  search  must 
be  sufficiently  fine  to  assure  that  the  iterative  procedure  does  not  converge  to  a  local  maximum. 


3.5.1  Sub-optimal  Technique 

The  coarse  parameter  search  and  iterative  procedure  described  above  can  be  computationally 
expensive  for  long  filter  lengths.  If  the  filter  is  of  length  P,  each  step  in  the  parameter  search  will  require 
the  estimation  and  inversion  of  a  P  by  P  covariance  matrix.  The  computational  cost  in  terms  of  complex 
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multiplications  is  on  the  order  of  F^.  Thus  anything  that  can  be  done  to  reduce  P  will  result  in  substantial 
computational  savings. 

One  method  to  reduce  the  effective  filter  size  is  weight  pruning  [66].  This  is  a  process  whereby 
only  the  largest  weights  in  the  filter  are  included  in  the  matrix  inversion  and  subsequent  computations. 
The  basic  idea  is  to  include  the  taps  which  correspond  to  strong  scattering  while  neglecting  the  taps 
corresponding  to  weak  scattering.  To  determine  which  taps  to  include,  one  would  apply  (3.20)  without 
any  Doppler  terms  to  a  very  short  sequence  of  multipath  data  and  retain  some  number  of  the  largest  taps. 
For  data  sequences  on  the  order  of  1-3  ms,  the  data  may  be  considered  stationary  since  the  differential 
Doppler  frequencies  are  relatively  small.  The  relative  amplitude  of  the  taps  obtained  in  this  manner  will 
be  close  to  the  amplitude  of  the  taps  in  the  full  configuration  when  the  correct  Doppler  terms  are  applied 
to  the  filter.  If  the  number  of  taps  can  be  reduced  from,  for  example,  100  to  10,  the  computation  savings 
will  be  on  the  order  of  10^,  which  is  considerable. 

3.6  CHAPTER  SUMMARY 

This  chapter  contains  the  design  of  the  ML  estimator  based  on  an  assumption  that  the  bistatic 
clutter  is  homogeneous.  The  structure  of  the  estimator  was  shown  to  be  a  time- varying  FIR  filter  with 
Doppler  modulators  to  match  the  Doppler  shifts  in  the  multipath  signals.  The  performance  of  this 
estimator  will  be  studied  in  Chapter  5.  The  optimum  solution  requires  a  dense  multidimensional  search 
over  emitter  parameter  space.  A  sub-optimum  technique  that  employs  weight  thinning  has  been 
suggested.  This  technique  retains  only  the  taps  corresponding  to  the  strongest  scattering  cells.  This 
concept  leads  naturally  to  the  estimator  based  on  an  inhomogeneous  clutter  model  that  will  be  described 
in  the  next  chapter. 
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4.  LOCALIZATION  IN  INHOMOGENEOUS  CLUTTER 


4.1  INTRODUCTION 

In  this  chapter  an  estimator  is  designed  that  is  based  on  the  assumption  that  the  bistatic  clutter  is 
inhomogeneous  and  dominated  by  a  number  of  discrete  scatterers.  This  is  a  two-part  estimation 
procedure.  The  receiver  first  measures  the  scattering  parameters  of  relative  time  delay,  D,  differential 
Doppler,  /,  azimuth,  0,  and  amplitude,  a,  for  some  or  all  the  dominant  discrete  scattering  points.  The 
scattering  parameter  estimates  are  then  used  to  estimate  the  emitter  parameters  through  a  set  of  nonlinear 
bistatic  localization  equations.  These  nonlinear  equations  are  linearized  using  a  truncated  Taylor  series 
expansion  that  retains  only  the  terms  less  than  second  order.  This  permits  the  use  of  linear  estimation 
techniques  such  as  best  linear  unbiased  estimator  (BLUE)  [67].  This  estimator  will  be  shown  to  be  the 
MLE  when  the  observation  noise  is  Gaussian.  As  with  the  estimator  in  the  previous  chapter,  application 
of  this  estimator  requires  a  coarse  search  followed  by  an  iterative  procedure. 


4.2  BASIC  MODEL  FOR  INHOMOGENEOUS  SCATTERING 


The  fundamentals  of  bistatic  scattering  were  introduced  in  Section  2.3  and  a  homogeneous  clutter 
model  was  described  in  Section  3.2.  We  consider  in  this  section  a  clutter  model  that  is  inhomogeneous. 
This  model  assumes  that  the  bistatic  clutter  is  dominated  by  a  number  of  discrete  scatterers,  and  that  the 
emitter  multipath  signals  are  primarily  due  to  these  scatterers.  The  received  power  for  a'  discrete  scatterer 
was  given  as 


(2.3) 


where  cris  the  RCS  of  the  point  scatterer.  Sources  for  discrete  scattering  can  be  man-made,  such  as  water 
towers,  buildings,  etc.,  or  natural  terrain  features,  such  as  mountains,  cliffs,  etc.  An  example  of  a  discrete 
scattering  scenario  is  shown  in  Figure  4.1. 

The  RCS  of  the  discrete  scatterers  varies  from  scatterer  to  scatterer  with  small  scatterers  being 
more  frequent  than  large  scatterers.  A  model  commonly  used  for  higher  frequencies  is  shown  in  Table 
4.1  [56].  The  densities  in  this  table  are  higher  than  those  observed  by  the  Mountaintop  system  at  WSMR, 
where  the  density  is  closer  to  0.026  40-dBm^  scatterers  per  square  kilometer. 
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Figure  4.1.  Range  and  Doppler  contours  for  a  moving  emitter.  A  reference  beam  pointed  at  the  emitter  directly 
samples  the  emitter  waveform,  and  multipath  beams  collect  attenuated,  delayed,  and  Doppler-shifted  replicas 
scattered  by  the  discrete  scatterers. 


TABLE  4.1 

Discrete  Clutter  Model  [56] 


Radar  cross  section 

Density,  per  square  kilometer 

60  dBm^ 

0.026 

50  dBm® 

0.26 

40  dBm® 

2.6 
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4.3  STATISTICAL  FORMULATION 


The  statistical  formulation  and  design  of  the  estimator  consists  of  two  parts.  The  first  part  consists 
of  estimating  the  parameters  of  relative  time  delay,  D,  differential  Doppler,  o,  azimuth,  (p,  and  amplitude, 
a,  for  a  given  number  of  dominant  scatterers.  This  is  described  in  Section  4.3.1  below.  The  second  part 
consists  of  using  the  scattering  parameters  to  estimate  the  emitter  parameters  of  range,  R,  heading,  H, 
velocity,  V,  and  altitude,  A.  This  will  be  described  in  Section  4.3.2. 

4J.1  Part  1:  Estimation  of  Scattering  Parameters 

Each  discrete  scatterer  that  reflects  the  emitter  signal  is  described  by  four  parameters  that  must  be 
estimated.  The  necessary  parameters  are  scatterer  delay,  Doppler,  azimuth,  and  amplitude  which,  for  the 
i"'  scatterer,  are  written  as  [Dj,  (pi,  a,]. 

An  estimator  based  on  cross-correlation  is  shown  in  Figure  4.2  for  a  single  scatterer.  This  structure 
is  the  maximum  likelihood  estimator  for  narrowband  signals  where  both  the  signal  and  the  additive  noise 
are  white.  We  search  over  the  independent  parameters  of  [D,  O),  for  the  maximum  cross-correlation 
IRszi.  The  joint  estimation  of  differential  time  delay  and  differential  Doppler  shift  was  shown  to  be 
uncorrelated  in  [43]  when  the  observation  time  is  sufficiently  long  and  the  signal-to-noise  ratios  are 
sufficiently  high.  Under  these  assumptions  the  joint  estimator  is  efficient  in  the  sense  that  it  is  unbiased 
and  achieves  the  Cramer-Rao  lower  bound.  We  assume  that  these  conditions  hold  here.  The  amplitude 
estimate,  a,  can  be  derived  from  the  cross-correlation,  r„,  as  a,.  =  {i)l  (0) . 

If  there  are  P  scatterers  in  the  data,  there  will  be  P  peaks  in  the  delay-Doppler-azimuth  space 
covered  by  the  estimator.  We  choose  the  parameters  at  these  peaks  to  be  our  scatterer  estimates.  To 
simplify  the  analysis  of  the  scatterer  estimates,  we  assume  that  each  scatter  can  be  treated  as  a  separate 
estimation  problem  with  no  interaction  between  neighboring  scatterers.  This  may  be  justified  somewhat 
by  the  fact  that  there  are  three  factors  working  to  decorrelate  the  scattered  signals:  delay,  Doppler,  and 
azimuth. 

With  the  above  assumptions,  the  variances  of  the  scatterer  parameter  estimates  for  this  estimator 
are  taken  directly  from  the  Cramer-Rao  lower  bounds  (CRLB)  for  point  targets.  We  are  assuming  that 
discrete  scatterers  have  the  same  bounds  as  conventional  radar  point  targets.  The  CRLB  for  delay, 
Doppler,  angle,  and  amplitude  will  be  given  in  Chapter  5. 
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Figure  4.2,  Maximum  likelihood  estimator  for  independent  scattering  parameters  [D,  (o,  0  ]. 


4,3.2  Part  2:  Estimation  of  Emitter  Parameters 


In  this  section  we  design  an  estimator  for  emitter  parameters  0=  [/?,  H,  V,  A]  based  on  the  scatterer 
parameter  vectors  [D,  that  were  determined  above.  We  begin  by  writing  the  measured  Doppler 

shift  for  the  /'*  scatterer  as  a  function  g  of  the  true  emitter  parameters  0and  the  remaining  true  scattering 
parameters  a,.  =  [D,  ,0,  ]  as 

S,=g(e,aJ+e^^  ,  (4.1) 


where  is  the  error  in  the  i""  Doppler  frequency  estimate.  The  function  g  is  a  highly  nonlinear  function 
that  is  described  in  Appendices  A,  B,  and  C  for  different  geometrical  approximations.  We  do  not  know 
the  true  scatterer  parameters  CCi,  so  we  expand  (4.1)  in  a  truncated  Taylor  series  about  the  estimates  as 


<0,  =  g{e,di )+ {d,d, )+ ie,d, )+ , 


where  and  (Z,-2  denote  the  delay  and  azimuth  estimates,  respectively.  We  are  assuming  that  the 

errors  in  scattering  parameters  are  small  and  have  retained  only  the  zero-  and  first-order  terms.  This 
linear  truncation  of  the  Taylor  series  permits  the  use  of  a  best  linear  unbiased  estimator  BLUE  [67] 
approach.  Note  that  -CC-d^  and  g'^^  denotes  the  partial  derivative  of  g  with  respect  to  the  variable 

Oil. 
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We  next  assume  that  we  have  an  initial  guess  that  is  close  to  the  true  emitter  parameters  which  are 
written  as  6  =  {61,62, 6^,0^}  =  {R,H,V,A} .  This  starting  point  is  denoted  as  ^  and  will  initially  be 

obtained  using  a  coarse  least-squares  search.  This  will  be  discussed  in  more  detail  in  Section  4.3.4  below. 
We  further  expand  (4.2)  in  a  Taylor  series  about  ft  as 

=  g'^  i^^,d.}+Sg  g'g  (4-3) 


2  2 


S  g  (0  ,d.)+5  g  is  ,6t]+e 

a  V  0  i'  0  * (o. 

i\  il  i2  t2  i , 

which  linearizes  the  problem  in  So.  The  term  in  square  brackets  denotes  the  noise  caused  by  errors  in 
estimating  the  scatterer  delay,  azimuth,  and  Doppler,  and  will  be  referred  to  as  observation  noise 

« .  =  si,  (^0  A)+  8  a,,  (^0  ’«/  )  +  ^co,  ■ 

To  determine  6e  we  rewrite  (4.3)  as 

which  can  be  written  in  vector  form  as 

AS =z-n  , 


(4.4) 


(4.5) 


(4.6) 


with 


8e2^o^^\) 

A  = 

8ei  802  (®o’^2 )  802  804  (^o’^2  ) 

•  •  •  • 

s= 

8e,(^oAp)  Se^i^oAp)  8e3(^oAp)§e4(^oAp)_ 

5.. 

and 


z  = 


'd)i-g(0o,dj' 

^ttii  (^0  ’^1  )■*"  ^021  ^a2i  (^0  ’^1  )“*"  ^©1 

^2-8(OoA) 

n  = 

^ai2  ^ ai2  (^0  ’^2  )“*"  ^^22  ^022  (^0  ’^2  )“*■  ^©2 

S^-g(6o,d^l 

^a:p8aip  ^a2p8a2p  )“*■  ^(Op  _ 

Equation  (4.6)  is  in  the  form  of  a  BLUE  (best  linear  unbiased  estimator),  which  is  solved  for  5  as 

S  =  [a^r;'a]-‘a^r;'z  , 

which  has  minimum  variance  of 

var(ft)  =  var(<^,.)  =  [(a^R^'a)  % 


(4.7) 


(4.8) 


45 


and  a  parameter  covariance  matrix  of 


R,  =  (a^r;'a)-‘ 


(4.9) 


The  covariance  matrix  of  the  observation  noise  Rn  is 

r„=r,+hXh  (4.10) 

with 


Rn 

0  R, 


and 


Sd 


H"  = 


0 


(4.11) 


(4.12) 


The  covariances  for  the  individual  scattering  estimates  are 

.2 


R.= 


0 


Rn  = 


.2 

D2 

0 


R.  = 


<  0 
0  •. 


(4.13) 


(4.14) 


(4.15) 


which  are  assumed  to  be  known  a  priori  for  this  development.  This  assumption  is  standard  for 
localization  problems  [22].  When  we  apply  this  technique  to  Mountaintop  data  in  Section  6.3,  we  will 
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use  the  CRLB  expressions  that  will  be  given  in  Chapter  5,  Equations  (5.14)-(5.16),  along  with  estimates 
of  SNR  for  each  scatterer,  to  determine  the  variance  of  the  scatterer  parameters. 

Note  that  the  PDF  of  the  measurement  errors  was  not  included  in  any  of  the  above  development, 
only  the  covariance.  This  is  a  property  of  the  BLUE  estimator.  When  the  PDF  of  the  measurement  errors 
is  Gaussian,  this  estimator  is  also  the  MVU  (minimum  variance  unbiased)  estimator  as  well  as  the  ML 
(maximum  likelihood)  estimator. 

4  J.3  Physical  Interpretation 

We  note  the  similarity  between  the  BLUE  estimator  described  by  (4.7)  and  the  weighted  least 
squares  estimator  (LSE)  [67]  which  can  be  written  as 

5  =  [a^Wa]''a^Wz.  (4.16) 

The  weighting  matrix  W  is  equal  to  the  inverse  of  measurement  error  covariance  matrix  R„  described  by 
(4.10) 

w  =  r;;'  .  (4.17) 

With  the  weighted  LSE  interpretation,  we  can  examine  (4.16)  in  more  detail  to  see  what  intuition  can  be 
gained  for  the  localization  problem  at  hand.  The  measurement  error  covariance  matrix  R„  is  the  sum  of 
the  covariance  matrix  for  the  Doppler  measurement  Ra,  and  the  product  H^RaH.  The  matrix  H  is  a  matrix 
of  gradients,  and  Ra  contains  the  covariance  matrices  of  the  time  delay  measurements  R^  and  the  angle 
measurements  R0. 

Consider  the  contribution  from  both  a  strong  scatterer  and  a  weak  scatterer.  The  high  SNR  for  the 
strong  scatterer  results  in  smaller  measurement  errors  than  that  of  the  weak  scatter.  Consequently,  the 
weighting  matrix  W  will  have  a  large  weight  corresponding  to  strong  scatterer  and  a  smaller  weight 
corresponding  to  the  weak  scatterer.  This  is  a  fundamental  property  of  weighted  LS  estimation. 

The  influence  of  the  two  matrices  of  gradients,  H  and  A,  is  not  immediately  obvious,  however. 
Consider  two  scatterers  with  equal  signal-to-noise  ratios  but  unequal  gradient  terms  in  H.  The  scatterer 
possessing  the  smaller  gradient  will  be  weighted  more  heavily  in  the  weighting  matrix  W.  The 
interpretation  is  that  the  smaller  gradient  means  that  a  given  measurement  error  will  produce  a  smaller 
change  in  the  desired  parameters  S.  The  opposite  is  true  for  the  gradient  matrix  A  that  describes  changes 
in  scatterer  differential  Doppler  as  a  function  of  the  desired  emitter  parameters.  If  two  scatterers  are 
otherwise  equal  but  one  of  them  possesses  a  larger  gradient  in  A,  then  the  contribution  of  that  scatterer 
will  be  weighted  more  heavily.  Such  a  scatterer  can  be  thought  of  as  possessing  more  information  about 
the  emitter  parameters.  This  effect  will  sometimes  allow  weaker,  distant  scatterers  to  count  more  heavily 
in  the  solution  than  strong  close-range  scatterers. 
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4.3.4  Solution  Technique 


The  estimation  procedure  begins  by  first  estimating  the  differential  delay,  Doppler,  and  azimuth  of 
P  dominant  scatterers  as  described  in  Section  4.3.1.  With  the  scatterer  parameters  available,  we  proceed 
with  the  development  in  Section  4.3.2,  which  requires  an  initial  guess  at  the  emitter  parameters  ft.  This 
can  be  obtained  by  performing  a  coarse  multidimensional  least-squares  search  over  the  parameter  space 
of  ft  which  can  be  written  as 

P  2 

=  •  (4.18) 

1=1 

As  with  the  estimator  developed  in  the  previous  chapter,  we  can  omit  emitter  altitude  from  the  coarse 
search  if  desired.  Any  nominal  altitude  will  be  close  enough  to  start  the  iterative  procedure  required  for 
the  final  estimates. 

We  denote  the  result  of  the  coarse  search  as  ft  and  determine  fusing  (4.7).  We  then  replace 

- ^0+^  (4.19) 

and  repeat  the  computation.  We  repeat  the  iterations  until  the  norm  of  8  is  sufficiently  small.  When  this 
occurs,  the  error  in  the  emitter  parameters  is  described  by  (4.8)  above. 

4.4  CHAPTER  SUMMARY 

A  two-part  estimator  has  been  designed  that  is  based  on  an  inhomogeneous  clutter  assumption  The 
scatterer  parameters  of  relative  time  delay,  differential  Doppler,  azimuth,  and  amplitude  were  estimated 
first.  These  scatterer  estimates  were  then  used  to  estimate  the  desired  emitter  parameters  of  range, 
heading,  velocity,  and  altitude.  The  estimator  was  shown  to  be  the  MVU  and  ML  estimator  when  the 
PDF  of  the  measurement  noise  is  Gaussian.  The  eovariance  of  the  emitter  parameters  was  also 
determined.  This  estimator  will  be  used  to  localize  emitters  using  simulations  in  Chapter  5  and  using 
Mountaintop  field  data  in  Chapter  6. 
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5.  ESTIMATOR  PERFORMANCE 


5.1  INTRODUCTION 

The  theoretical  performance  of  the  estimators  designed  in  Chapters  3  and  4  is  considered  in  this 
chapter  via  the  Cramer-Rao  lower  bounds  (CRLB).  The  CRLB  describe  the  minimum  variance  for 
parameter  estimates  that  an  unbiased  estimator  can  achieve.  The  CRLB  will  be  determined  for  each 
estimator  and  evaluated  for  a  number  of  localization  scenarios.  The  effect  of  varying  the  emitter 
parameters  of  range,  heading,  velocity,  altitude,  and  power,  as  well  as  the  receiver  parameters  of 
bandwidth,  beamwidth,  steering  direction,  integration  time,  etc.,  will  be  studied.  We  also  consider  the 
consequences  of  some  of  the  simplifying  assumptions  that  were  used  to  design  the  estimators  and  how 
these  assumptions  could  lead  to  optimistic  performance  calculations  when  applied  to  experimental  field 
data.  The  estimator  based  on  the  homogeneous  clutter  iriodel  is  discussed  in  Section  5.2  and  the  estimator 
based  on  the  inhomogeneous  clutter  model  is  discussed  in  Section  5.3.  Monte-Carlo  simulations  are 
performed  for  each  estimator  to  validate  both  the  bounds  calculations  and  the  overall  implementation  of 
the  estimators. 

5^  ESTIMATOR  FOR  THE  HOMOGENEOUS  CLUTTER  MODEL 


5.2.1  Determination  of  the  Cramer-Rao  Lower  Bounds 


The  covariance  of  the  emitter  paranieter  estimates  is  equal  to  the  inverse  of  the  Fisher  information 
matrix  [68], [69].  The  diagonal  elements  of  the  inverted  matrix  are  equal  to  the  CRLB.  The  elements  of 
the  Fisher  information  matrix  are  defined  by 


LV 


de. 


A 


dO-dOj 


\nf^(z\e)\  , 


(5.1) 


where  is  the  likelihood  function  that  was  developed  in  Chapter  3  for  the  homogeneous  clutter 

model.  The  log  likelihood  function  is  written  as 


ln/^(z|0)=ln 


f  1 


—  N—2N 

K  <J 


■((z-m(e)f  {z-mie))/a^)  . 


Using  the  right  side  of  (5.1)  we  write 

le 


F,=-E\ 


^In (z|0)|  =  ~ >n(0))" (z - m(0))]/<7' |  . 


(5.2) 


(5.3) 


The  product  in  square  brackets  can  be  written  as 
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z  +|m(0)|  -m(0)^z--z^m(0)  , 


which  can  be  expanded  in  terms  of  the  real  and  imaginary  parts  of  m  and  z  as 
X  \^r  (")  +  zf  (n)  +  in)  +  mf  (n)  -  2(m^  {n)z,  (n)  +  m,.  (n)z,.  («)))  . 

w=l 

We  insert  (5.5)  into  (5.3)  to  obtain  the  elements  of  the  Fisher  information  matrix, 


F,..  =  ^E 


X  {^r  +  zf  +  +  m?  -  2(m^Zr  +  miZi)) 


1  rJ  f  V  ^  ^  „  dm: 

=  2m,  — +2m,^-2z,35^-2z,^ 


Mj  -  36, 


J  V  d^m^  d^m^  _  d^m^  d\ 

MS  dOj  de^ddj  ddidOj  deidej~^‘  dOiddj 


_  2  ^  dnij.  dm^  dm^ 

The  last  simplification  made  use  of  the  fact  that  E(z}  =  m.  (The  sample  index  n  has  been  omitted  from 
each  term  for  convenience  but  we  note  that  the  sum  is  over  N,  the  entire  length  of  the  data  vector.)  This 
result  is  similar  to  the  results  for  frequency  estimation  described  elsewhere  [70].  Equation  (5.6)  can  also 
be  written  in  vector  form 

2  dm{6)"  dm(d) 

^  de,  ■ 


To  make  use  of  this  result  in  (5.6),  we  rewrite 
m{n)=Y,  Op  s{n  - 


P 

min)  =  in)  +  jm.  («)  =  X  IK”  “  'C>p  )|  cos(d)^  («  -  I>p )  +  (p,^„_D  )+<?«,) 

p=i  p  " 

P  (5.8) 

+7Xklk(«-^p)|sin(®p(n-£>p)  +  <?),(„_o  )  +<?>„,)  > 

p=\  P 

which  uses  an  amplitude-phase  representation  for  the  complex  filter  coefficients  ap  and  emitter  waveform 
s(n).  To  compute  the  Fisher  information  matrix  we  use  (5.8)  in  (5.6)  and  differentiate  with  respect  to  the 
desired  parameters.  (Note  that  the  emitter  waveform  s(n)  is  included  in  the  bound  calculation.)  The  three 
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parameters  we  are  mainly  interested  in  are  emitter  range  R,  heading  H,  and  velocity  V.  However,  we  see 
from  (5.8)  that  each  emitter  multipath  beam  will  contain  2P  nuisance  parameters  corresponding  to  the 
amplitude  and  phase  of  each  filter  coefficient.  Thus  the  Fisher  information  matrix  will  be  a  square  matrix 
with  2BP+3  rows  and  columns  where  B  is  the  number  of  beams.  The  emitter  parameters  of  (R,  H,  V)  are 
contained  in  the  (Op  term  which  is  given  by 

/  \  iTtV  (  .  -\(  (Op  +  ^ 

©  {r,H,v)  = -  cosff-cos  H  +  sm  — r — v,  Ti  ^  ‘  ^ 

Wp\i^,n,y )  ^  dI+2[R^ +RDp}^-cos^)\\ 

y  \  y  /  / J 

for  the  two-dimensional  approximation  that  ignores  altitude.  This  equation  is  derived  in  Appendix  A. 

The  details  of  the  calculation  of  the  Fisher  information  matrix  are  relegated  to  Appendix  D; 

however,  the  form  of  the  matrix  is  interesting  and  worth  a  few  comments.  If  we  write  ¥■■  —  > 

then,  after  we  perform  the  differentiation  and  expectation  operations,  we  are  left  with  a  matnx  (for  a 
single  beam)  of  the  form 
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(5.10) 


for  the  case  of  the  filter  delays  spaced  equal  to  or  greater  than  the  correlation  length  of  the  emitter 
waveform  s(n).  The  zeros  in  the  matrix  indicate  a  lack  of  coupling  between  parameter  estimates,  i.e.. 
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uncertainty  in  one  parameter  estimate  does  not  affect  the  variance  that  can  be  achieved  for  the  other 
parameter  estimate.  For  this  problem,  we  can  see  that  the  amplitude  estimates  are  uncoupled  from  each 
other  and  from  all  other  parameters.  However,  the  phase  estimates  are  coupled  with  the  desired 
parameters  and  the  desired  parameters  are  coupled  with  each  other.  This  means,  for  example,  that  our 
ability  to  measure  emitter  range  when  the  emitter  velocity  is  unknown  will  be  degraded  from  the  case 
when  the  emitter  velocity  is  known. 


Due  to  the  block  structure  of  the  Fisher  information  matrix,  it  will  not  be  necessary  to  invert  the 
whole  2P +3  matrix.  We  can  write 


(5.11) 


where  B  is  the  F  by  F  diagonal  matrix  of  derivatives  with  respect  to  the  scatterer  amplitudes.  Hence  the 
inverse  of  the  information  matrix  is 


F-‘ 


A'  0 

0 


(5.12) 


and  we  need  only  to  evaluate  A  *.  For  the  variance  of  the  emitter  parameters,  the  only  terms  of  interest 
are  the  first  three  diagonal  elements  of  A’\ 

If  we  want  to  estimate  the  emitter  altitude  as  well  as  range,  velocity,  and  heading,  then  the  size  of 
t  e  information  matnx  will  be  2FF+4.  A  three-dimensional  approximation  for  the  Doppler  terms  that 
includes  the  emitter  altitude  is 
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where  A  is  the  emitter  altitude  and  F„  is  the  range  from  the  emitter  to  the  scattering  patch.  This  equation 
IS  denved  in  Appendix  B.  It  reduces  to  (5.9)  when  the  emitter  altitude  is  set  to  zero. 

In  the  remainder  of  the  section  we  shall  compute  the  CRLB  for  several  engagement  scenarios.  We 

shall  also  perform  Monte-Carlo  simulations  which  serve  to  validate  the  bounds  calculations  and 
simulations. 


5.2,2  CRLB  Calculations  and  Behavior 

The  performance  of  these  estimators  is  scenario  dependent.  The  CRLB  depend  on  both  the  emitter 
parameters  and  the  receiver  parameters.  In  this  section  we  study  the  behavior  of  the  bounds  while  varying 
several  of  these  parameters.  The  emitter  parameters  that  we  shall  consider  include  power,  range,  heading. 
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velocity,  and  altitude.  For  the  receiver,  we  consider  bandwidth,  beamwidth,  steering  direction,  numba  of 
beams,  integration  time,  etc.  While  this  section  focuses  on  the  estimator  for  the  homogeneous  clutter  model, 
most  of  these  results  are  directly  applicable  to  the  estimator  for  the  inhomogeneous  clutter  model  that  is 
covered  in  Section  5.3. 

Once  the  localization  scenario  and  filter  parameters  are  specified,  the  CRLB  are  determined  by 
inverting  (5.10).  The  initial  localization  scenario  is  shown  in  Figure  5.1  with  the  corresponding  delay  and 
Doppler  contours  spaced  by  25  |lis  and  50  Hz,  respectively.  The  emitter  range  is  100  km,  the  altitude  is  7 
km,  and  the  velocity  is  150  m/s  with  a  heading  of  90  degrees.  The  nominal  receiver  azimuth  is  30  degrees 
and  the  other  receiver  parameters  were  taken  from  Table  2.1.  The  nominal  emitter  power  is  500  Watts 
spread  over  the  200  kHz  bandwidth  of  the  receiva. 


Figure  5.1.  Basic  localization  scenario  for  CRLB  analysis. 


The  tapped-delay-line  architecture  of  the  localization  filter  was  shown  in  Figure  3.3.  For  the  bounds 
calculations,  the  filter  employs  50  coefficients  (or  taps)  spaced  5  |xs  apart  for  a  total  delay  extent  of  250  ms 
(or  equivalently,  75  km).  The  magnitude  of  the  filta:  coefficients  is  determined  using  the  constant-gamma 
scattering  model  of  Equation  (3.2)  with  a  gamma  of  -20  dB.  Figure  5.2(a)  shows  a  map  of  the 
bistatic  reflectivity  in  the  8-degree  receive  beam.  The  average  value  is  about  -45  dB  meters  per 
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square  meter.  Figure  5.2(b)  shows  the  receive  power  per  square  meter  of  terrain.  The  corresponding 
emitter  signal-to-noise  ratio,  versus  delay  and  Doppler,  is  shown  in  Figure  5.2(c)  for  an  integration  period 
of  40  ms.  This  interval  corresponds  to  16  PRIs  of  recorded  Mountaintop  data.  The  amplitude  of  the  filter 
coefficients  is  shown  in  Figure  5.2(d).  This  is  the  ratio  between  the  scattered  power  and  the  direct-path 
power  and  is  hence  independent  of  the  emitter  transmit  or  receiver  noise  power. 


Figure  5.2.  Outputs  from  scattering  model  include  (a)  reflectivity  map,  (b)  scattered  power  map,  (c)  delay-Doppler 
map,  and  (d)  filter  coefficient  amplitude  for  an  8-degree  beam  at  30  degrees  azimuth. 


54 


The  first  parameter  we  shall  consider  is  the  emitter  velocity.  Figure  5.3(a)-(d)  illustrates  the  RMS 
error  (or  standard  deviation)  calculations  for  the  four  localization  parameters  of  range,  heading,  velocity, 
and  altitude  as  a  function  of  emitter  velocity.  We  are  showing  the  square  root  of  the  emitter  variance 
given  by  the  Cramer  Rao  lower  bounds.  We  see  from  the  figure  that  the  emitter  range,  heading,  and 
altitude  estimates  are  inversely  proportional  to  emitter  velocity.  When  the  emitter  velocity  is  zero,  the 
range,  heading,  and  velocity  errors  are  infinite  and  this  localization  technique  fails.  (Hence  the  inclusion 
of  the  word  “moving”  in  the  report  title.)  Note  that  the  RMS  velocity  errors  are  independent  of  the 
velocity  of  the  emitter.  A  similar  result  occurs  in  tone  estimation  problems,  i.e.,  the  ability  to  estimate  the 
frequency  of  a  tone  is  generally  independent  of  the  tone’s  frequency. 

We  next  compute  the  CRLB  as  a  function  of  emitter  altitude.  The  results  are  shown  in  Figure  5.4. 
The  Doppler  expressions  were  taken  from  the  three-dimensional  approximation  described  in  Appendix  B. 
We  see  that  the  ability  to  estimate  the  emitter  parameters  of  range,  heading,  and  velocity  is  essentially 
independent  of  emitter  altitude,  but  that  the  RMS  altitude  error  is  inversely  proportional  to  the  emitter 
altitude.  This  behavior  can  be  understood  in  terms  of  conical  angles  versus  azimuthal  angles.  As  the 
emitter  altitude  increases,  the  conical  angles  depart  more  quickly  from  the  azimuthal  angles  and  this 
facilitates  better  altitude  estimation.  Another  point  is  that  the  RMS  errors  for  the  other  parameters  can  be 
seen  to  increase  slightly  as  the  altitude  increases.  This  behavior  is  explained  by  the  fact  that  as  the 
altitude  increases,  the  Doppler  shift  at  each  patch  on  the  ground  decreases  slightly.  The  decrease  in 
Doppler  is  equivalent  to  a  decrease  in  emitter  velocity  which,  from  Figure  5.3,  was  shown  to  increase 
parameter  errors. 

The  CRLB  versus  receiver  azimuth  is  shown  in  Figure  5.5  and  the  CRLB  versus  emitter  heading  is 
shown  in  Figure  5.6.  The  behavior  of  the  bounds  in  these  figures  is  more  difficult  to  interpret  without 
looking  at  the  individual  expressions  for  the  elements  in  the  Fisher  information  matrix.  One 
generalization  that  can  be  made  is  that  the  more  Doppler  spread  across  the  length  of  the  filter,  the  better 
the  estimator  will  perform. 

There  are  several  additional  parameters  upon  which  the  localization  performance  will  depend. 
These  include  the  scattered  SNR,  integration  time,  number  of  receive  beams,  receiver  bandwidth, 
beamwidth,  etc.  Some  of  the  relations  may  be  observed  by  examining  the  expressions  for  the  CRLB  or  be 
obtained  by  evaluating  the  bounds  directly.  Many  of  these  additional  parameters  will  be  discussed  in  the 
remainder  of  this  subsection.  Some  of  the  relations  cannot  be  determined  fully  from  the  CRLB  described 
here  due  to  assumptions  and  simplifications  made  in  the  statistical  formulation  of  the  likelihood  function. 
Examples  of  this  are  the  effect  of  changing  the  receive  bandwidth  or  beamwidth  and  these  will  be 
described  below  as  well. 

The  variance  of  the  emitter  parameter  estimates  is  inversely  proportional  to  SNR  or,  alternatively, 
the  RMS  error  is  inversely  proportional  to  the  square  root  of  the  SNR.  The  SNR  of  the  scattered  signals 
is  proportional  to  the  product  of  emitter  power  and  the  terrain  reflectivity  as  well  as  the  area  of  the  delay- 
azimuth  cell.  The  SNR  is  (approximately)  inversely  proportional  to  the  square  of  the  emitter  range. 
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Figure  5.4.  CRLB  for  the  emitter  parameters  of  (a)  range,  (b)  heading,  (c)  velocity,  and  (d)  altitude  all  as  a 
function  of  emitter  altitude. 
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The  variance  of  the  parameter  estimates  is  inversely  proportional  to  the  sample  length  (or 
integration  time)  to  the  third  power.  This  is  similar  to  tone  estimation  problems  [70]  where  the  variance 
of  the  frequency  estimates  is  inversely  proportional  to  N^.  However,  as  integration  times  become  very 
lengthy,  Doppler  dispersion  may  be  an  important  factor.  That  is,  the  assumption  that  all  of  the  scattering 
in  a  single  delay-azimuth  cell  will  occur  at  the  same  Doppler  frequency  may  not  be  valid.  This  could  be 
circumvented  by  expanding  the  localization  filter  to  include  multiple  Doppler  modulators  per  delay  tap. 

Estimation  performance  improves  when  more  beams  are  included.  If  we  use  two  very  close  beams 
at,  say  30  and  30.1  degrees,  the  parameter  variance  is  reduced  by  a  factor  of  two.  If,  however,  we  use  two 
widely  spaced  beams  at  +30  and  -30  degrees,  the  variance  was  found  to  drop  by  a  factor  of  100.  For  the 
first  case,  very  little  new  information  is  provided  by  the  second  beam  at  30.1  degrees.  The  factor-of-two 
improvement  is  due  to  an  increased  gain  against  noise.  However,  for  the  second  case,  the  -30  degree 
beam  provides  a  different  look  at  the  emitter,  and  the  two  beams  combined  provide  a  significant 
performance  enhancement  over  either  single  beam  individually.  This  behavior  is  related  to  GDOP  which 
was  described  in  Section  2.4.2. 

We  next  consider  the  bandwidth  of  the  receiver.  The  reciprocal  of  the  receiver  bandwidth  was  used 
to  specify  tap  spacing  in  the  localization  filter.  However,  if  we  reduced  the  bandwidth  and  corresponding 
number  of  taps  by  a  factor  of  two,  the  performance  remains  unchanged  because  the  size  of  the 
delay-azimuth  cell  (and  corresponding  scattered  SNR)  increases  by  a  factor  of  two.  The  net  effect,  based 
on  the  statistical  formulation  presented  herein,  is  no  change  in  performance.  There  is  a  hidden  cost  to  a 
reduction  in  bandwidth  however.  As  the  area  of  the  delay-azimuth  cell  increases,  the  amount  of  Doppler 
dispersion  in  a  given  cell  will  increase.  As  was  the  case  with  increased  integration  times,  multiple 
Doppler  modulators  may  be  required  to  cope  with  the  increase  in  Doppler  dispersion. 

The  same  effect  occurs  for  an  increase  in  the  receive  beamwidth.  If  the  receive  antenna  beamwidth 
increases,  its  gain  will  drop  but  the  product  of  gain  and  the  beamwidth  is  constant.  Thus  there  will  be  no 
change  in  the  CRLB  as  calculated  above,  but  again  the  increased  cell  area  means  increased  Doppler 
dispersion.  The  effect  of  Doppler  dispersion  on  the  bounds  will  be  discussed  further  in  Section  5.4. 

We  next  present  the  results  of  some  Monte-Carlo  simulations  and  follow  with  a  discussion  of  the 
performance  of  the  estimator  that  is  based  on  an  inhomogeneous  clutter  model.  That  estimator  has  more 
concrete  ties  to  receiver  bandwidth,  beamwidth,  and  integration  time. 

5.2.3  Monte-Carlo  Simulations 

Monte-Carlo  simulations  were  performed  to  assure  that  the  estimates  generated  by  the  ML 
estimator  described  in  Chapter  3  achieve  the  CRLB  determined  in  this  chapter.  The  simulations  consisted 
of  generating  data  vectors  using  Equation  (3.3)  for  many  different  realizations  of  the  emitter  waveform 
and  noise,  applying  the  localization  filter  of  Figure  3.3  while  keeping  track  of  the  statistics  of  the 
estimates.  This  provides  a  good  check  on  the  derivation  and  implementation  of  the  model,  estimator,  and 
CRLB  calculations.  However,  the  simulations  do  not  necessarily  describe  how  well  the  model  will  match 
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field  Hata  collected  by  an  actual  receiver  over  real  terrain.  Application  of  the  CRLB  to  real  data  will  be 
considered  further  in  Section  5.4. 

Before  proceeding,  we  first  review  the  properties  of  maximum  likelihood  estimators  (taken  from 

[17]): 

1.  The  MLE  is  asymptotically  unbiased.  This  means  that  as  the  number  of  observations  increase,  the 
bias  tends  to  zero. 

2.  The  MLE  is  consistent.  This  means  that  the  mean-squared  estimation  error  tends  to  zero  as  the 
number  of  observations  increases. 

3.  The  MLE  is  asymptotically  efficient.  An  efficient  estimator  is  one  that  achieves  the  CRLB. 

4.  Asymptotically,  the  maximum  likelihood  estimate  is,  in  general,  distributed  as  a  Gaussian  random 
variable. 

This  section  examines  properties  1,  3,  and  4.  (Property  2  was  addressed  in  the  previous  section  where  we 
noted  that  the  CRLB  was  inversely  proportional  to  the  third  power  of  the  number  of  observations.)  These 
properties  state  that  if  we  perform  enough  independent  simulations,  we  should  obtain  parameter 
estimates  that  are  Gaussian  distributed  (property  4)  with  means  equal  to  the  true  parameters  (property  1) 
and  variances  equal  to  the  CRLB  (property  3). 

The  radar-emitter  scenario  for  the  Monte-Carlo  simulations  is  the  same  as  that  used  in  Section 
5.2.2.  The  CRLB  versus  emitter-to-noise  ratio  was  computed,  and  over  100  simulations  were  performed 
for  each  of  four  different  emitter  power  levels.  The  results  are  shown  in  Figure  5.7.  The  asterisks 
represent  the  RMS  estimation  error  for  the  independent  simulations.  We  see  from  the  figure  that  the 
estimator  achieves  the  CRLB  as  expected.  (Some  statistical  fluctuation  around  the  CRLB  can  also  be 
expected.) 

Another  way  to  view  the  simulation  results  is  to  plot  the  density  functions  of  the  parameter 
estimates.  As  noted  above,  the  estimates  should  be  Gaussian  random  variables  with  variances  described 
by  the  CRLB.  Figure  5.8  shows  the  density  of  the  emitter  range  and  heading  estimates  as  histograms 
along  with  normalized  Gaussian  PDFs  shown  for  reference.  The  normalized  Gaussian  PDFs  use  the  true 
parameters  of  range  or  heading  for  the  mean  and  the  CRLB  for  the  variance.  The  two  figures  show  that 
Monte-Carlo  estimates  fit  the  Gaussian  distributions  well. 
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Figure  5.7.  Comparison  of  Monte-Carlo  simulations  (*)  to  CRLB  calculations  for  (a)  range  and  (b)  heading. 
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Figure  5.8.  (a)  Emitter  range  and  (b)  heading  histograms  for  110  simulations.  The  mean  range  was  99.85  and  the 
mean  heading  was  90.16.  The  solid  curve  represents  the  theoretical  PDF  based  on  the  CRLB  calculations 


5.3  ESTIMATOR  FOR  THE  INHOMOGENEOUS  CLUTTER  MODEL 


The  covariance  of  the  emitter  parameters  for  the  estimator  based  in  an  inhomogeneous  clutter 
model  was  specified  in  Chapter  4  by  Equations  (4.9)-(4.15).  These  expressions  relate  the  localization 
errors  to  the  errors  present  in  the  scatterer  delay,  Doppler,  and  azimuth  estimates  through  the  gradient 
matrices  A  and  H.  Since  these  matrices  were  already  specified  in  Chapter  4,  all  we  need  specify  to 
complete  the  bounds  calculations  is  the  variance  of  the  scattering  parameters. 

The  variance  of  the  scattering  parameters  will  be  taken  from  the  CRLB  for  point  targets.  Hence  we 
are  assuming  that  the  discrete  scatterers  present  in  the  inhomogeneous  clutter  have  the  same  bounds  as 
conventional  point  targets.  The  bounds  depend  on  SNR  as  well  as  characteristics  of  the  receiver 
including  bandwidth,  beamwidth,  and  integration  time.  The  CRLB  for  delay  D,  Doppler  /,  angle  0,  and 
amplitude  a  can  be  shown  to  be,  respectively  [67], [71] 
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where  J]  =  is  the  mainbeam  SNR  per  sample,  B  is  the  signal  bandwidth,  T  is  the  signal 

duration,  M  is  the  number  of  elements  in  the  array,  and  <(>  is  the  angle  to  the  scatterer  measured  from 
boresight.  The  product  t}N  represents  the  SNR  at  the  output  of  the  correlator  which  has  been  increased  by 
a  factor  of  N  over  the  per-sample  value.  Note  that  the  SNR  T]  is  actually  an  estimate  f]  = 

because  it  contains  the  amplitude  estimate  d .  The  variance  of  the  additive  thermal  noise  is  assumed 

to  be  known.  Typical  values  for  the  Mountaintop  system  for  a  point  scatterer  with  an  SNR  of  10  dB 
are  Ob=0.26  km,  G(^2\A  rad/s,  and  cr^O.Tl  deg.  In  general,  the  emitter  localization  errors  for  the 
Mountaintop  system  will  be  dominated  by  the  Doppler  errors  which  are  larger  in  magnitude  than  the 
delay  or  angle  terms. 

5.3.1  Parameter  Covariance  Calculations 

The  localization  parameter  covariances  can  be  determined  once  the  scatterer  SNR  and  scatterer 
spatial  distribution  have  been  specified.  For  a  general  simulation,  we  would  select  scatterers  from  a 
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uniform  distribution  as  suggested  in  Section  2.2.1.  However,  we  begin  here  by  selecting  a  scatterer  SNR 
and  distribution  that  matches  the  scattering  levels  for  the  homogeneous  case  discussed  in  Section  5.2.2. 
That  is,  we  locate  50  scatterers  on  the  30-degree  receive  beam,  evenly  spaced  in  bistatic  delay,  with  SNR 
levels  corresponding  to  the  levels  present  for  the  homogeneous  case.  This  facilitates  a  direct  comparison 
between  the  parameter  variances  for  the  inhomogeneous  case  and  the  parameter  variances  for  the 
homogeneous  case  when  the  scattering  levels  are  assumed  to  be  equivalent. 

As  expected,  the  result  for  the  two  calculations  are  essentially  identical  as  can  be  seen  in  Figure 
5.9.  Although  both  estimators  use  the  data  differently,  we  achieve  the  same  result  because  both  ML 
estimators  minimize  the  variance  of  the  emitter  parameters.  The  first  estimator  chooses  the  set  of  emitter 
parameters  that  deterministically  produces  the  set  of  Doppler  modulators  that  best  match  the  received 
data,  while  the  second  estimator  selects  emitter  parameters  that  best  match  the  set  of  discrete  scatterer 
parameters.  The  likelihood  function  for  the  homogeneous  case  deals  directly  with  the  measured  data 
while  the  likelihood  function  for  the  inhomogeneous  case  is  in  terms  of  the  scatterer  estimates. 


The  small  differences  apparent  in  Figure  5.9  are  probably  due  to  the  different  approximations 
made  during  the  statistical  formulation  for  each  case.  For  example,  the  first  estimator  treated  time  delay 
and  azimuth  as  known  independent  variables  while  the  second  estimator  treated  these  as  noisy  estimates. 
Conversely,  the  second  estimator  treats  the  scattering  SNR  as  if  it  were  known  exactly,  while  the  first 
estimator  treats  scattering  amplitudes  as  nuisance  parameters.  The  likelihood  function  for  the 
inhomogeneous  estimator  can  be  written  in  terms  of  the  Doppler  measurement  vector  c5  and  the  emitter 
parameters  0as 
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for  the  simplified  case  where  the  time-delay  and  azimuth  errors  can  be  ignored. 

The  fact  that  the  two  estimators  can  be  shown  to  be  identical  when  the  scattering  is  the  same 
underscores  the  fact  that  most  of  the  observations  that  were  made  in  Section  5.2.2  also  apply  for  this 
estimator.  That  is,  this  estimator  will  show  the  same  dependence  on  emitter  velocity,  SNR,  etc.,  as  the 
estimator  for  the  homogeneous  clutter  model  did. 


For  the  second  example,  we  use  the  same  scatterer  SNR  and  delay  as  the  previous  example,  but 
select  the  azimuth  of  the  scatterers  from  a  uniform  random  distribution  between  ±45  degrees.  The  result 
is  shown  in  Figure  5.10  and  exhibits  considerable  improvement  over  the  previous  case.  This  performance 
improvement  is  also  related  to  GDOP.  That  is,  if  we  consider  each  scatterer  to  be  a  virtual  sensor,  then 
better  performance  is  achieved  by  spreading  the  sensors  over  the  terrain  than  by  combining  them  to  a  line 
at  a  particular  azimuthal  angle. 
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Figure  5.10.  CRLB  versus  velocity  for  50  scatterers  spead  in  azimuth  from  -45  to  +45. 


5.3.2  Monte-Carlo  Simulation 


A  Monte-Carlo  simulation  was  performed  for  the  previous  example  where  the  scatterer  azimuth 
was  selected  from  a  uniform  distribution.  The  scatterer  parameters  of  delay,  Doppler,  and  azimuth  were 
corrupted  from  their  true  values  by  additive  noise  with  variance  equal  to  the  CRLB  expressions  of 
(5.14)-(5.16).  The  results  are  shown  in  Figure  5.11.  As  expected,  the  estimates  are  unbiased  and  match 
the  theoretical  variances. 
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Velocity  (m/s) 


Figure  5.11.  Histograms  for  1700  simulations.  The  mean  of  the  simulations  equal  99.991  km,  90.000  deg., 
180.005  m/s,  7.015  km  for  range,  heading,  velocity,  and  altitude,  respectively.  The  solid  curve  represents  the 
theoretical  PDF  based  on  the  CRLB  calculations. 
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5.4  APPLICATION  OF  CRLB  TO  REAL  DATA 


Sections  5.2  and  5.3  contained  theoretical  bounds  calculations  based  on  homogeneous  and 
inhomogeneous  clutter  assumptions.  An  important  question  is  whether  or  not  these  bounds  apply  to 
parameter  estimates  obtained  using  experimental  field  data.  The  answer  depends  on  the  capabilities  of 
the  receiver  and  the  scattering  properties  of  the  terrain.  We  now  consider  the  two  estimators  separately. 

The  statistical  formulation  used  for  the  estimator  for  homogenous  clutter  assumed  narrow  receive 
beams  and  closely  spaced  delay  cells.  This  permitted  the  use  of  a  single  Doppler  term  to  describe  the 
scattered  signal  from  each  delay-azimuth  cell.  For  nonzero  beamwidth  and  delay  cells,  some  Doppler 
dispersion  will  be  present  and  this  will  reduce  the  performance  of  the  estimator  from  the  theoretical 
results.  Another  factor  that  could  reduce  performance  of  this  estimator  is  the  clutter  itself.  If  the  clutter  is 
not  homogeneous,  but  instead  contains  large  clutter  discretes,  emitter  energy  could  enter  the  filter  via  the 
beam  edges  and  sidelobes.  The  effect  of  this  energy  would  be  to  increase  the  apparent  noise  floor  and 
thereby  reduce  the  performance  of  the  estimator. 

The  key  assumption  for  the  inhomogeneous  clutter  estimator  is  that  the  clutter  behaves  as  a 
collection  of  independent  discrete  scatterers.  We  assumed  that  we  could  estimate  the  scatterer  parameters 
as  well  as  we  can  estimate  point  radar  targets.  This  model  will  be  violated  somewhat  by  clutter  that  is 
more  homogeneous  in  nature.  In  that  case,  the  dominant  scatterers  could  possess  more  spatial  extent  and 
be  spread  in  any  combination  of  delay,  azimuth,  or  Doppler.  The  parameter  estimates  for  extended 
scatterers  will  have  larger  errors  th^ul  point  scatterers  with  equivalent  SNR. 

The  factors  described  above  can  be  thought  of  as  modeling  error.  That  is,  our  models  do  not 
describe  the  physical  problem  completely.  The  signal  components  resulting  from  modeling  error  is 
considered  model  noise.  Thus  the  total  noise  g  is  the  sum  of  the  receiver  thermal  noise  gr  and  the  model 
noise  g„.  This  will  be  discussed  further  in  Chapter  6  where  we  shall  see  that  the  model  noise  for  the 
Mountaintop  system  parameters  is  10-20  dB  for  the  homogeneous  model  and  less  than  10  dB  for  the 
inhomogeneous  model.  The  difference  in  the  levels  of  model  noise  for  these  two  models  suggests  that 
the  WSMR  terrain  is  more  inhomogeneous  than  homogeneous.  Nevertheless,  we  shall  see  in  the  next 
chapter  that  both  estimators  will  successfully  localize  real  emitters  contained  in  Mountaintop  data 
collected  at  WSMR. 

5.5  CHAPTER  SUMMARY 

The  theoretical  performance  of  the  two  estimators  designed  in  Chapters  3  and  4  was  determined  in 
this  chapter.  The  behavior  of  the  bounds  was  considered  for  a  number  of  localization  scenarios.  The  two 
estimators  were  demonstrated  to  be  equivalent  for  identical  scattering  scenarios.  Monte-Carlo 
simulations  were  provided  to  substantiate  the  localization  and  bounds  calculations.  In  the  next  chapter  we 
will  utilize  these  techniques  with  Mountaintop  field  data. 
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6.  APPLICATION  TO  MOUNTAINTOP  DATA 


6.1  INTRODUCTION 

In  this  chapter  we  use  the  estimators  developed  thus  far  to  localize  several  airborne  emitters 
contained  in  Mountaintop  field  data.  This  data  was  collected  at  White  Sands  Missile  Range  (WSMR) 
where  the  terrain  is  a  mixture  of  deseat  and  mountains.  Both  the  estimator  for  homogeneous  clutter  and  the 
estimator  for  inhomogeneous  clutter  will  be  used.  Estimation  errors  obtained  with  real  data  will  be 
compared  to  the  CRLB  when  appropriate.  This  chapter  will  serve  to  demonstrate  the  techniques  presented 
in  the  previous  three  chapters. 

The  Mountaintop  system  was  described  briefly  in  Chapter  2.  A  photograph  of  the  antenna  and  the 
airborne  emitter  is  shown  in  Figure  6.1.  For  the  localization  work  presented  ha*e,  the  antenna  was  rotated 
90  degrees  to  the  horizontal  polarization  mode  referred  to  as  RSTER-90.  In  this  mode,  the  anteima 
supports  adaptability  in  the  azimuthal  dimension.  Thus,  once  the  data  for  each  element  have  been  recorded, 
receive  beams  can  be  formed  in  any  or  aU  directions  simultaneously.  The  emitter  contains  dipole  antennas 


Figure  6.1.  Photograph  of  Mountaintop  installation  and  Lear  jet  emitter  at  North  Oscura  Peak,  NM. 
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in  the  nose  and  tail  portions  of  the  airplane.  The  radiation  patterns  for  these  dipole  antennas  are  fairly 
broad.  The  nose  antenna  is  used  when  the  emitter  is  flying  inbound  toward  the  receiver  and  the  tail  is 
used  when  flying  outbound  away  from  the  receiver.  The  location  of  the  receiver  and  the  emitter  for 
several  data  sets  is  shown  in  Figure  6.2.  The  emitter’s  ranges  vary  from  78.3  km  for  HOT-6161  to 
212  km  for  HOT-6132.  Additional  information  on  the  Mountaintop  system  and  experiments  can  be  found 
in  [6], [7], [9]. 


Figure  6.2.  Map  showing  Mountaintop  radar  site  and  location  of  emitters. 


Three  emitters  will  be  localized  using  the  two  estimators.  The  data  sets  are  referred  to  as  HOT- 
6067,  HOT-6161,  and  HOT-6132  and  were  collected  on  April  6,  1994  for  the  purposes  of  studying 
terrain  scattered  interference  [72].  The  estimator  for  homogeneous  clutter  will  be  employed  in  Section 
6.2  and  the  estimator  for  inhomogeneous  clutter  will  be  employed  in  Section  6.3.  In  Section  6.4  we 
consider  the  case  of  reduced  emitter  power  and  bandwidth.  We  summarize  the  results  of  this  chapter  in 
Section  6.5. 
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6.2  RESULTS  USING  THE  ESTIMATOR  FOR  HOMOGENEOUS  CLUTTER 


6.2.1  Data  Set  HOT-6067 

The  first  real  emitter  to  be  localized  using  out-of-plane  multipath  was  that  of  data  set  HOT-6067. 
This  emitter  was  located  101.7  km  west  of  RSTER  traveling  northeast  at  a  velocity  of  184.5  m/s.  The 
engagement  geometry  is  shown  in  Figures  6.2  and  6.3.  The  true  emitter  parameters  were  obtained  from  a 
GPS  receiver  located  on  board  the  Lear  jet.  The  emitter  azimuth  was  determined  by  selecting  the  peak  of 
a  conventional  (nonadaptive)  beamformer  and  found  to  be  within  about  0.1  degrees  of  that  indicated  by 
the  GPS.  The  coordinate  system  used  to  describe  emitter  position  in  this  chapter  is  shown  in  Figure  6.4. 
This  coordinate  system  measures  angles  relative  to  the  receiver  boresight.  The  full  three-dimensional 
coordinate  system,  including  emitter  altitude,  is  shown  in  Appendix  C. 


deg.  @  1 84.5  m/s 


Figure  6.3.  Geometry  for  data  set  HOT-6067.  Emitter  heading  and  Mountaintop  boresight  are  given  in  degrees 
relative  to  true  north. 


Range  R 


Heading  H 
Velocity  V 


Figure  6.4.  Emitter  parameter  estimates  of  heading  and  velocity  are  determined  relative  to  the  receive  antenna 
boresight. 

The  general  estimator  described  in  Chapter  3  specifies  that  multiple  receive  beams  be  used  to 
localize  an  emitter.  However,  for  this  case,  a  single  beam  was  all  that  was  required.  (This  is  a  function  of 
the  emitter  parameters,  in  particular,  the  emitter  heading.)  The  localization  filter  used  75  taps  spaced  by  4 
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fis  with  data  selected  from  a  75  ms  coherent  processing  interval  (CPI).  The  receive  beam  was  pointed  to  - 
42  degrees  which  is  out  in  front  of  the  emitter.  A  four-dimensional  grid  search  was  performed  over  the 
emitter  parameters  of  range,  heading,  velocity,  and  altitude.  Prior  to  performing  the  grid  search,  the  data 
are  first  passed  through  a  3-tap  version  of  the  filter  to  cancel  the  emitter  signal  entering  the  receive  beam 
via  the  antenna  siddobes.  Strictly  speaking,  only  one  tap  should  be  required  for  such  a  sidelobe  canceUCT. 
The  additional  taps  are  used  to  cancel  the  strong  multipath  signals  that  originate  from  the  glistening  surface 
and  from  scattera-s  very  close  to  the  receiver. 

The  output  of  the  four-dimension  search  is  the  four-dimensional  log-likelihood  function  (with 
coordinates  of  range,  heading,  velocity,  and  altitude)  that  we  desire  to  maximize.  To  visualize  this  surface, 
we  show  several  cuts  or  slices  through  it.  The  range-heading-velocity  cut  is  shown  in  Figure  6.5.  For  this 
cut,  the  altitude  estimate  was  fixed  at  7.7  km  which  corresponds  to  the  peak  of  the  Mkelihood  function.  We 
see  a  single  dominant  peak  near  the  true  emitter  parametors  and  some  smallCT  peaks  and  other  structure 
throughout.  A  range-altitude  slice  though  the  likdihood  fimction  is  shown  in  Figure  6.6.  The  emitter 
heading  and  velocity  were  fixed  at  the  values  corresponding  to  the  likelihood  peak.  The  one-dimensional 
cuts  through  the  peak  are  shown  in  Figure  6.7  (a)-(d),  respectively,  for  range,  heading,  velocity,  and 
altitude.  Tlie  dashed  lines  denote  the  true  emitter  parameters. 


Figure  6.5.  3-D  log-likelihood  function  for  HOT-6067  at  a  fixed  altitude. 
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EMITTER  RANGE  (km) 


Figure  6.6.  Range-altitude  slice  for  fixed  heading  and  velocity. 


The  localization  results  are  listed  in  Table  6.1  along  with  the  true  parameter  values  and  the 
estimation  orors.  The  earors  are  listed  by  their  actual  value  and  in  terms  of  the  CRLB  drived  in  Chapter 
5.  One  can  argue  that  if  the  estimation  OTors  are  less  than  4-5  o,  then  we  are  operating  close  to  the  CRLB. 
Hence,  for  this  example,  we  have  far  exceeded  the  CRLB  for  range  and  heading  estimates.  Howeva:,  we 
had  pointed  out  in  Section  5.4  that  the  CRLB  was  developed  based  on  the  assumption  that  the  data  match 
the  model  perfectly.  That  is,  the  clutter  is  homogenous  (without  clutter  discretes)  and  described  by  a 
constant  gamma  clutter  model.  Different  values  of  gamma  will  produce  different  CRLB.  We  also  assumed 
that  the  enargy  scattered  by  each  delay-azimuth  cell  is  modulated  by  a  single  Doppler  shift  term  The 
presence  of  clutter  discretes  in  the  beam  edges,  or  Doppler  dispersion  due  to  the  extended  nature  of  the 
scattering  cells,  will  produce  model  noise.  We  will  see  in  Section  6.4  that  it  is  possible  to  determine  the 
amount  of  irKxlel  noise  present  by  obsaving  the  residual  noise  level  at  the  output  of  the  localization  filter 
when  it  is  operating  with  the  correct  emitter  parameters.  For  this  data  set,  the  model  noise  was  found  to  be 
20  dB  above  the  receiver  thermal  noise.  A  20-dB  increase  in  noise  increases  the  CRLB  RMS  error 
calculations  by  10  dB  and  thus  all  the  estimates  are  within  3a  of  the  true  values. 


73 


(a) 


(b) 


(d) 


Figure  6.7.  Slices  through  the  four-dimensional  log-likelihood  function  at  the  emitter  (a)  range,  (b)  heading,  (c) 
velocity,  and  (d)  altitude.  The  dashed  lines  denote  the  true  parameter  values  obtained  via  GPS. 
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TABLE  6.1 

Estimation  Results  for  Data  Set  HOT-6067 


Emitter 

Parameter 

True 

Value 

Estimated 

Value 

Error 

Value  CRLB 

CRLB 

a 

Range  (km) 

101.7 

97 

4.7 

22a 

0.21 

Heading  (deg) 

-147.2 

-149.5 

2.3 

29a 

0.078 

Velocity  (m/s) 

184.5 

185.5 

1.1 

2.2a 

0.5 

Altitude  (km) 

8 

7.7 

.3 

2.7a 

0.11 

6.2.2  Data  Set  HOT-6161 

The  next  data  set  we  will  examine  is  HOT-6161.  The  position  of  the  emitter  was  shown  in  Figure 
6.2.  For  this  data  set,  the  emitter  was  flying  very  low  through  an  area  known  as  Milligan  Gulch.  The 
point  of  the  Milligan  Gulch  experiment  was  to  try  to  deny  the  receiver  the  direct-path  line-of-sight  signal 
in  order  to  stress  TSI  mitigation  algorithms  that  rely  on  the  direct  path  signal.  The  emitter  altitude  varied 
from  500  to  1000  feet  above  the  ground  but  the  direct-path  signal  was  not  severely  attenuated  (probably 
due  to  diffraction). 

For  this  data  set,  we  use  several  beams  to  perform  the  localization.  The  emitter  power  versus 
azimuth  is  shown  in  Figure  6.8.  The  conventional  nonadaptive  beamformer  and  MVDR  (minimum 
variance  distortionless  response)  [4], [73]  beamformer  are  both  used.  The  conventional  beamformer  is 
useful  for  measuring  the  SNR  of  the  emitter  direct  path  signal  but  fails  to  measure  the  emitter  multipath 
power  at  other  angles  due  to  the  high  sidelobes.  Conversely,  the  MVDR  spectrum  captures  the  multipath 
signal  versus  angle  but  fails  to  measure  the  power  in  the  emitter  direct  path.  The  later  point  is  probably 
due  to  small  system  calibration  errors  or  the  correlation  in  the  multipath  signals  or  some  combination  of 
both  of  these  effects. 

Five  receive  beams  were  formed  to  collect  emitter  multipath  energy.  The  beam  azimuths  are  given 
by  the  vertical  lines  in  the  figure.  These  beams  were  selected  because  they  correspond  to  higher  than 
average  emitter  multipath.  Cuts  through  the  likelihood  functions  for  each  beam  are  shown  in  Figure 
6.9(a)  and  (b)  for  altitude  and  range,  respectively.  The  altitude  cuts  have  been  normalized  for  unity  peak 
in  the  likelihood  function  while  the  relative  magnitude  of  the  range  cuts  has  been  preserved.  The  largest 
peak  is  from  the  62-degree  beam  while  the  shortest  peak  is  from  the  20-degree  beam.  Part  (b)  of  the 
figure  shows  that  all  beams  produce  peaks  near  the  correct  emitter  range.  In  part  (a)  we  see  that  all  beams 
produce  a  peak  (albeit  broad)  near  the  true  emitter  altitude. 
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Figure  6.8.  Emitter  SNR  versus  receiver  azimuth  for  HOT-6161  using  conventional  and  MVDR  beamformers. 
Emitter  and  beam  azimuthal  angles  are  shown  as  vertical  lines. 


The  final  estimates  are  taken  by  forming  the  sum  (or  mean)  of  the  individual  likelihood  functions  for 
the  five  beams.  The  emitta  range-heading  slice  is  shown  in  Figure  6.10  and  the  estimated  quantities  and 
estimation  errors  are  listed  in  Table  6.2. 


TABLE  6.2 

Estimation  Results  for  Data  Set  HOT-6161 


Emitter 

Parameter 

True 

Value 

Estimated 

Value 

Error 

Value  % 

Range  (km) 

78.3 

83.0 

mm 

6% 

Heading  (deg) 

121.8 

122.1 

0.25% 

Velocity  (m/s) 

134.3 

139 

D 

3.6% 

Altitude  (km) 

0.5 

0 

- 
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ALTITUDE  (km) 

(a) 


RANGE  (km) 

(b) 


Figure  6.9.  (a)  Altitude  and  (b)  range  cuts  through  the  log-likelihood  functions  for  five  multiple  beams.  (Data 
HOT-6161.) 
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Figure  6.10.  Range-heading  slice  through  four-dimensional  log-likelihood  function  for  HOT-6161. 


6.2.3  Data  Set  HOT-6132 


Several  additional  data  sets  were  localized  with  results  similar  to  the  previous  two  sections.  The 
third  data  set  discussed  in  this  chapter  is  HOT-6132.  This  set  is  of  particular  interest  because  the 
estimator  successfully  used  above  fails  to  localize  this  emitter.  The  range  for  this  emitter  is  212  km  and 
the  geometry  was  shown  in  Figure  6.2. 

Ten  evenly  spaced  receive  beams  were  formed  covering  from  -60  to  +60  degrees  of  azimuth.  The 
range  cuts  for  the  ten  beams  are  shown  in  Figure  6.1 1.  We  can  see  four  or  five  peaks  close  to  the  true 
emitter  range  but  no  consensus  between  the  results  and  hence  this  localization  technique  has  failed  in  this 
case. 


The  principal  cause  for  this  failure  is  evident  in  the  delay-Doppler  map  shown  in  Figure  6.12.  The 
solid  line  represents  the  delay-Doppler  ridge  for  the  -60  degree  beam  and  the  dashed  lines  denote  the  Vi 
power  points.  For  this  case,  there  is  clearly  a  great  deal  of  Doppler  spread  across  the  beam  at  any  given 
delay.  For  example,  at  200  p,s,  the  spread  covers  over  100  Hz  while  the  Doppler  resolution  is  only  6.5- 
19.5  Hz  depending  upon  the  length  of  the  coherent  processing  interval  (CPI).  Multipath  signals 
originating  off  the  ridge  by  more  than  one  half  the  reciprocal  of  the  integration  time  will  essentially  be 
decorrelated  from  the  Doppler  modulator  in  the  filter  and  will  not  produce  a  peak  in  the  likelihood 
function.  Instead,  these  signals  manifest  themselves  as  increases  in  the  noise  level.  Also  evident  in  the 
figure  is  the  general  inhomogeneous  nature  of  the  bistatic  clutter  for  this  beam.  While  some  of  the  energy 
is  confined  to  the  delay-Doppler  ridge,  there  is  some  strong  scattering  outside  of  the  ridge  that  is  picked 
up  by  the  beam  edges.  (The  null-to-null  beamwidth  of  this  beam  is  close  to  40  degrees.)  There  is  also  a 
very  strong  scatterer  located  at  15  jxs  that  produces  significant  delay  sidelobes.  All  of  these  factors  reduce 
the  effectiveness  of  this  filter  to  the  point  where  it  did  not  produce  the  correct  estimate. 

There  are  several  possible  ways  one  could  try  to  improve  this  situation.  One  approach  would  be  to 
build  an  antenna  with  a  narrower  beamwidth.  If,  for  example,  a  1-2  degree  wide  beam  could  be  formed, 
this  narrow  beam  would  act  to  confine  the  energy  to  the  Doppler  ridge  and  reduce  Doppler  dispersion. 
Another  approach  to  capture  the  Doppler  spread  is  to  add  several  Doppler  modulators  per  delay  cell.  This 
would  increase  the  dimensions  of  the  filter  and  thereby  increase  the  computational  burden.  Another 
possible  approach  could  be  to  reduce  the  integration  time  and  thereby  decrease  the  Doppler  resolution. 
However,  we  saw  in  Chapter  5  that  localization  performance  is  proportional  to  the  cube  of  the  integration 
time. 


Instead  of  trying  to  fix  this  estimator  that  was  designed  to  work  in  homogeneous  clutter,  we  now 
switch  to  the  estimator  designed  for  inhomogeneous  clutter.  The  clutter  inhomogeneities  that  caused 
problems  for  the  estimator  used  in  this  section  will  promote  more  accurate  results  for  the  estimator  used 
in  the  next  section. 
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RELATIVE  DOPPLER  (Hz) 


Figure  6.12.  Delay-Doppler  map  shows  wide  Doppler  spread  across  (-60  deg)  beam. 
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6.3  RESULTS  USING  THE  ESTIMATOR  FOR  INHOMOGENEOUS  CLUTTER 

In  this  section  we  demonstrate  the  estimator  that  was  designed  for  inhomogeneous  cluttea-  using  the 
same  three  data  sets  examined  in  the  previous  section.  We  shall  see  that  this  estimator  performs  better  than 
the  estimator  described  above.  This  is  probably  due  to  the  general  inhomogeneous  nature  of  the  bistatic 
clutter  at  WSMR. 

The  estimator  designed  for  inhomogeneous  clutter  is  a  two-part  procedure  where  first  the  scatta-er 
parameters  are  estimated  and  then  the  emitter  parameters  are  estimated.  The  basic  ML  estimator  for 
scatterer  parameters  was  shown  in  Figure  4.2.  This  structure  requires  that  we  search  for  scatterers  over  all 
angles,  delays,  and  Doppler  bins.  The  practical  implementation  of  this  scatterer  parameter  estimator  is 
shown  in  Figure  6.13.  Cross-correlations  were  performed  for  the  direct  path  and  thirty-three  beams  ranging 
firom  -64  to  64  degrees.  In  this  manner,  a  three-dimensional  cube  of  scattared  energy  versus  delay, 
Doppler,  and  a/imnth  is  formed.  These  correlation  cubes  have  been  generated  in  the  past  by  Gabel  [66]. 
The  cube  is  passed  to  the  “Select  and  Centroid  Scattarers”  box  which  applies  a  threshold  and  centroids  the 
results  to  a  set  of  discrete  scatterer  parameter  vectors  D,  (»,  (|),  and  S/N.  The  threshold  is  set 
significantly  high  so  that  false  alarms  are  essentially  nonexistent. 


Figure  6.13.  Estimator  for  scatterer  parameters  of  delay,  Doppler,  azimuth,  and  SNR. 


Some  of  the  problems  the  “Select  and  Centroid  Scatterars”  box  has  to  contend  with  are  depicted  in 
Figure  6.14.  This  figure  shows  a  delay-Doppler  map  for  a  typical  multipath  beam.  In  the  map  one  can  see 
sidelobes  in  all  three  dimensions  of  delay,  Doppler,  and  azimuth.  This  box  rejaits  sidelobes  by  selecting 
only  the  largest  scatterer  for  each  delay,  Doppler,  and  azimuth  dimensioa 
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Figure  6.14.  Delay-Doppler  map  showing  energy  in  the  delay,  Doppler,  and  azimuth  sidelobes.  The  delay 
sidelobes  are  the  horizontal  bands  of  energy  and  the  Doppler  sidelobes  are  vertical  bands.  Scatterers  outside  of 
the  dashed  lines  are  off  of  the  delay-Doppler  ridge  and  are  therefore  in  the  beam’s  edges  or  sidelobes. 


Another  potential  problem  is  that  estimates  for  closely  spaced  discrete  scatterors  that  are  not 
separated  by  several  resolution  cells  in  delay,  Doppler,  or  azimuth,  may  be  biased  due  to  coherent 
interference.  This  is  dqpicted  graphically  in  Figure  6.15  and  is  independent  of  SNR.  Scatterer  parameter 
errors  of  this  nature  are  not  accounted  for  in  the  scatterer  CRLB  computations  and  therefore  the  bounds 
calculations  will  be  optimistic.  Howevo",  as  we  shall  see  below,  this  estimator  produces  fairly  accurate 
parameter  estimates  that  are  close  to  the  CRLB  in  several  cases. 
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(a)  Scatterers  in  phase 


Figure  6.15.  Examples  of  errors  produced  by  closely  spaced  scatterers.  True  system  responses  and  measured 
responses  for  (a)  in-phase  and  (b)  out-of-phase  combinations  are  shown.  The  two  in-phase  scatterers  merge  into  a 
single  extended  scatterer  while  the  out-of-phase  scatterers  seem  to  “repel”  each  other. 


6.3.1  Data  Set  HOT-6067 

The  scatterer  estimates  for  HOT-6067  are  shown  in  Figure  6.16.  The  three-dimensional  cube  is 
collapsed  to  azimuth-delay  in  part  (a)  and  azimuth-Doppler  in  part  (b).  The  location  of  each  scatterer  or 
scattering  center  is  denoted  by  an  ‘x’.  Short-delay  scatterers  (<  10  |js)  and  scatterers  within  a  beamwidth  or 
so  of  the  emitter  azimuth  are  not  used.  These  are  believed  to  be  biased  by  the  strong  direct  path  signal.  The 
total  number  of  scatterers  used  for  the  localization  of  this  data  set  is  67. 
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Figure  6.16.  Cross-correlation  plots  for  HOT-6067.  Part  (a)  shows  correlation  vj.  time,  part  (b)  shows  correlation 
vs.  Doppler,  and  part  (c)  shows  correlation  vs.  delay  and  Doppler  for  a  single  beam  at  -40  degrees.  The  position 
of  each  scatterer  is  denoted  by  an  ” 

With  the  scatterer  parameters  estimated,  we  are  now  ready  to  estimate  the  emitter  parameters.  A 
coarse  least-squares  search  using  Equation  (4.18)  is  performed  ova:  emitter  range,  heading,  and  velocity 
for  a  nominal  fixed  altitude.  The  range-heading  slice  is  shown  in  Figure  6.17.  The  coarse  search  is 
performed  to  provide  an  initial  guess  at  the  emitter  parameters.  This  initial  guess  is  used  as  the  starting 
point  for  the  iterative  procedure  described  in  Section  4.3.4.  We  found  that  it  was  not  necessary  to  include 
emitter  altitude  in  the  coarse  search.  The  iterative  procedure  converged  to  the  same  point  regardless  of  the 
initial  starting  altitude. 


40  60  80  100  120  140  160  180  200 
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Figure  6.17.  Least-squares  search  through  parameter  space  for  HOT-6067. 


The  iterative  procedure  convCTged  rapidly  to  the  estimates  listed  in  Table  6.3.  Only  five  iteration 
steps  were  required.  The  results  are  excellent  and  are  fairly  close  to  the  CRLB  for  all  parameters  except  the 
emitter  velocity.  The  range  error  is  only  0.7  km  and  the  heading  error  is  only  0.2  degrees. 


TABLE  6.3 

Results  for  Data  Set  HOT-6067 


Emitter 

Parameter 

True 

Value 

Estimated 

Value 

Error 

Value  CRLB 

CRLB 

c 

Range  (km) 

101.7 

101.0 

-0.7 

2.3c7 

0.31 

Heading  (deg) 

-147.2 

-147.0 

0.2 

4a 

0.05 

Velocity  (m/s) 

184.5 

179.3 

5.2 

11.1CT 

0.47 

Altitude  (km) 

8 

9.4 

1.4 

4.5ct 

0.31 
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6.3.2  Data  Set  HOT-6161 


The  same  approach  was  used  for  data  set  HOT-6161  with  the  same  excellent  results.  The  number 
of  scatterers  used  in  the  localization  was  61  and  the  localization  results  are  listed  in  Table  6.4.  As  with 
the  previous  case,  these  estimates  are  reasonably  close  to  the  CRLB.  This  fact  suggests  the 
inhomogeneous  model  contains  less  model  noise  and/or  estimation  biases  than  the  homogeneous  model 
when  applied  to  WSMR  terrain.  Model  noise  will  be  discussed  further  in  Section  6.4. 


TABLE  6.4 

Results  for  Data  Set  HOT-6161 


Emitter 

Parameter 

True 

Value 

Estimated 

Value 

Error 

Value  CRLB 

CRLB 

a 

Range  (km) 

78.3 

78.1 

0.2 

0.8o 

0.24 

Heading  (deg) 

121.8 

121.9 

3.3o 

0.03 

Velocity  (m/s) 

134.3 

134.9 

0.6 

3.5o 

0.17 

Altitude  (km) 

0.5 

-0.24 

1.6a 

0.46 

6.3.3  Data  Set  HOT.6132 

Data  set  HOT-6132  is  of  particular  interest  because  the  estimator  for  homogeneous  clutter  failed  to 
localize  the  emitter.  One  of  the  reasons  this  estimator  failed  is  the  inhomogeneous  nature  of  the  clutter  in 
this  data  set.  Hence  we  expect  the  performance  of  the  estimator  of  inhomogeneous  clutter  to  exceed  that 
of  the  estimator  for  homogeneous  clutter  for  this  case. 

Eighty  scatterers  were  used  for  the  localization  of  this  data  set.  The  results  of  the  coarse  search  are 
shown  in  Figure  6.18.  A  single,  well-defined  peak  is  present  in  the  figure.  Hence  this  estimator  has 
successfully  localized  this  emitter  and  the  results  are  listed  in  Table  6.5.  The  CRLB  for  this  case  are 
small  because  of  the  large  number  of  strong  scatterers  used  for  the  localization.  The  estimation  errors  in 
terms  of  CRLB  are  quite  large  compared  to  the  previous  two  examples.  With  the  exception  of  altitude, 
the  localization  errors  themselves  are  not  very  large,  however.  For  example,  the  range  error  is  only  3% 
but  the  18.4-km  (60,000  feet)  altitude  estimate  is  off  by  more  than  10  km. 
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Figure  6.18.  Least-squares  search  through  range-heading  parameter  space  for  HOT-6132. 


TABLE  6.5 

Results  for  Data  Set  HOT-'6132 


Emitter 

Parameter 

True 

Value 

Estimated 

Value 

Error 

Value  CRLB 

CRLB 

a 

Range  (km) 

211.9 

218.6 

6.7 

29o 

0.23 

Heading  (deg) 

-36.2 

-36.7 

-.5 

70o 

0.007 

Velocity  (m/s) 

142.1 

150.5 

8.4 

42a 

0.2 

Altitude  (km) 

8 

18.4 

10.4 

42a 

0.25 

In  an  attenpt  to  reduce  the  estimation  errors  for  this  data  set,  the  localization  was  repeated  for  one 
fewer  scatter.  The  localization  was  performed  80  times  while  the  contribution  of  each  of  the  80  scatterers 
was  omitted  one  at  a  time.  The  set  of  79  scatterers  that  resulted  in  the  lowest  altitude  estimate  was  selected 
as  the  new  estimate.  The  idea  was  that  perhaps  one  of  the  scatterer  estimates  was  severely  biased 
by  the  coherent  interference  effect  that  was  shown  in  Figure  6.15,  and  that  omitting  the  scatterers,  one  at 
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^  would  isolutc  the  bud  scuttercr.  The  new  estimutes  ure  shown  in  Xuble  6.6,  The  range  estimate  is 
now  within  the  CRLB  and  the  remaining  parameter  estimates  have  improved  as  well.  The  scatterer  that 
was  removed  had  an  SNR  of  44  dB  which  was  the  tenth  largest.  The  errors  are  still  larger  than  for  the 
previous  two  data  sets,  however. 

Another  possible  explanation  for  the  increased  errors  in  this  example  could  be  that  the  flat-earth 
assumption  is  being  stressed  for  this  longer  range  case.  Also,  perhaps  more  than  one  scatterer  estimate  is 
biased  by  the  coherent  interference  effect. 


TABLE  6.6 

Results  for  Data  Set  HOT-6132  With  One  Scatterer  Removed 


Emitter 

Parameter 

True 

Value 

Estimated 

Value 

Error 

Value  CRLB 

CRLB 

a 

Range  (km) 

211.9 

212.3 

0.4 

1.7a 

0.23 

Heading  (deg) 

-36.2 

-36.6 

-.4 

56c 

0.007 

Velocity  (m/s) 

142.1 

145.6 

3.5 

17.5a 

0.2 

Altitude  (km) 

8 

12.0 

4 

16a 

0.25 

6.4  RESULTS  FOR  REDUCED  EMITTER  POWER  AND  BANDWIDTH 

We  have  demonstrated  thus  far  the  ability  to  localize  a  strong  700  W  emitter’  in  a  fairly  strong 
clutter  environment.  What  about  low-power  emitters  or  weaker  clutter?  The  results  presented  in  Chapter 
5  suggest  that  it  should  be  possible  to  localize  emitters  with  much  lower  signal-to-noise  ratios.  To 
investigate  this  further  we  have  computed  the  theoretical  performance  limits  as  a  function  of  both  emitter 
power  and  clutter  strength.  We  considered  a  geometry  similar  to  the  HOT-6067  data  set  and  assumed  that 
all  the  energy  scattered  from  the  shaded  portion  of  Figure  6.19  could  be  collected  and  processed 
optimally.  The  result  is  shown  as  two  lines  in  Figure  6.20.  These  are  loci  of  constant  performance  as  a 
function  of  emitter  SNR  and  average  scattering  strength. 

Mountaintop  parameters  are  indicated  as  a  point  in  Figure  6.20.  This  is  the  theoretical  point  which 
assumes  that  the  models  are  perfect.  In  practice,  the  model  error  will  slide  this  point  to  the  right 
somewhat.  Nevertheless,  there  is  plenty  of  room  for  the  localization  of  much  weaker  emitters. 


'  The  emitter  power  is  2  kW  over  600  kHz  but  only  700  W  is  in  the  receiver  passband  of  200  kHz. 
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Figure  6.19.  Delay-Doppler  contours  for  the  multiple-beam  bounds  calculation.  The  area  that  is  shaded  is 
included  in  the  CRLB  calculation. 


Figure  6.20.  Loci  of  constant  performance  versus 


emitter  SNR  and  average  bistatic  clutter  strength. 


We  can  investigate  localization  performance  versus  SNR  with  Mountaintop  data  by  artificially 
increasing  the  noise  floor  of  the  receiver.  This  is  accomplished  by  the  addition  of  independent  noise 
sources  at  the  terminals  of  the  receiver  as  shown  in  Figure  6.21.  This  is  equivalent  to  reducing  either  the 
power  of  the  emitter  or  the  reflectivity  of  the  terrain.  Lowpass  filters  have  also  been  added  to  investigate 
localization  performance  as  a  function  of  emitter  bandwidth. 


Figure  6,21.  Emitter  localization  filter  for  showing  injected  noise  and  lowpass  filtering  operations. 


Emitter  localization  was  performed  for  HOT-6067  with  noise  added  in  10  dB  increments  from  0  to 
50  dB.  This  corresponds  to  emitter  signal-to-noise  ratios  of  60-10  dB  per  antenna  element  or  71.5-21.5 
dB  in  the  main  beam.  The  range  cuts  that  result  are  shown  in  Figure  6.22.  The  ordinate  axis  has  been 
changed  from  log-likelihood  to  cancellation  ratio  (Izl/lgl)  to  facilitate  an  estimation  of  the  model  noise 
which  is  discussed  below.  In  the  figure  we  see  peaks  near  the  correct  range  for  emitter  SNR  extending  all 
the  way  down  to  20  dB  per  element.  The  full  range-heading  slice  is  shown  in  Figure  6.23  for  both  the 
full-power  case  and  the  30-dB  SNR  case.  The  30-dB  SNR  case  corresponds  to  an  emitter  power  of  only 
0.350  Watts. 

We  can  infer  the  level  of  model  noise  for  this  scenario  by  considering  the  TSI  cancellation 
properties  of  the  localization  filter  and  the  results  shown  in  Figure  6.22.  We  see  in  the  figure  that  the 
10-dB  increase  in  receiver  noise  that  reduced  the  SNR  from  60  to  50  dB  produced  a  negligible  change  in 
cancellation  ratio.  This  is  a  clear  indication  that  the  cancellation  performance  of  this  filter  is  not  limited 
by  receiver  noise  but  by  model  noise.  The  second  10-dB  increase  in  noise  dropped  the  cancellation 
performance  by  about  3  dB  from  the  starting  point.  The  3-dB  drop  in  cancellation  indicates  that  this  is  the 
point  where  the  model  noise  and  the  receiver  noise  are  equal  in  power.  Hence  the  model  noise  for  this 
scenario  is  about  20  dB  above  the  receiver  noise.  This  suggests  that  the  theoretical  Mountaintop  point  on 
Figure  6.20  should  be  moved  to  the  left  by  20  dB  for  the  terrain  at  WSMR. 
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Figure  6.22.  Cancellation  ratio  of  localization  filter  versus  range  cut  for  emitter  SNR  ranging  from  60  down  to 
10  dB. 


The  effect  of  reducing  the  emitter  bandwidth  and  power  is  shown  in  Figure  6.23.  Part  (a)  shows  the 
full-power,  fuU-bandwiddi  case,  and  part  (b)  shows  the  30  dB  per  element  SNR  case  filtered  to  10-kHz 
bandwidth.  The  10-kHz  bandwidth  was  chosen  because  it  is  commensurate  with  the  bandwidth  of  a  voice 
communications  channd.  The  peak  of  the  likelihood  function  is  spread  somewhat  but  it  clearly  identifies 
the  emitter  very  near  the  true  range  and  heading.  The  number  and  spacing  of  the  localization  filter 
coefficients  was  not  changed  for  the  reduced  bandwidth  case. 

The  paformance  of  the  estimator  for  inhomogeneous  clutter  was  also  studied  as  a  function  of  emitter 
power  and  bandwidth.  The  results  are  listed  in  Table  6.7.  The  first  column  shows  the  results  presented 
above  for  the  full  powa*  and  full  bandwidth  case.  The  middle  column  shows  the  results  for  a  reduction  in 
power  of  43  dB.  The  number  of  scatterers  available  has  dropped  from  67  to  36  but  the  localization  results 
are  still  reasonable.  For  the  final  column,  the  emitter  power  was  dropped  by  33  dB  and  the  bandwidth 
reduced  to  20  kHz.  A  three-dimensional  estimate  was  made  that  may  still  be  useful.  For  this  final  case,  the 
number  of  scatterers  available  was  reduced  to  six.  FurthM'  reductions  in  SNR  or  bandwidth  caused  this 
estimator  to  faU  completely.  Thus  it  appears  that  the  estimator  for  homogenous  clutter  may  be  better  for 
lower  bandwidth  signals  than  the  estimator  for  inhomogeneous  clutter  for  this  data  set. 
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Figure  6.23.  Range-velocity  slices  of  the  log-likelihood  function  for  (a)  700W/200  kHz  and  (b)  0.35W/10  kHz 
cases. 


TABLE  6.7 

Results  for  Reduced  Emitter  Power  and  Bandwidth  Using  the  Estimator  for 

Inhomogeneous  Clutter 


700  W  Power 
200  kHz  BW 

0.035  W  Power 
200  kHz  BW 

0.35  W 
20  kh 

Power 

zBW 

Emitter 

True 

Est. 

CRLB 

Est. 

CRLB 

Est. 

CRLB 

Parameter 

Value 

a 

a 

a 

Range  (km) 

101.7 

101.0 

0.31 

94.0 

2.1 

78.8 

14 

Heading  (deg) 

-147.2 

-147.0 

0.05 

-146 

0.46 

-141 

4.3 

Velocity  (m/s) 

184.5 

179.3 

0.47 

175 

3.4 

130 

26.8 

Altitude  (km) 

8 

9.4 

0.31 

10.9 

0.35 

j 

— 

— 

#  Scatterers 

67 

36 

6 

6.5  CHAPTER  SUMMARY 

This  chapter  contains  a  demonstration  of  the  localization  techniques  designed  in  this  report  using 
Mountaintop  field  data.  The  localization  estimators  wae  found  to  perform  very  well  with  performance 
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reasonably  close  to  the  theoretical  performance  predicted  in  Chapter  5.  Localization  performance  for 
lower  power  and  bandwidth  emitters  was  also  considered  by  artificially  reducing  the  emitter  power  and 
bandwidth.  We  are  very  pleased  with  these  results  and  feel  that  we  have  fully  validated  the  estimators 
designed  in  Chapters  3  and  4.  We  consider  some  basic  extensions  to  these  results  in  the  next  chapter. 
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7.  EXTENSIONS  TO  THE  BASIC  RESULTS 

7.1  INTRODUCTION 

In  this  chapter  we  indicate  how  the  basic  results  of  Chapters  3-6  may  be  extended  to  include  a 
broader  range  of  localization  scenarios  and  class  of  emitter  signals.  Examples  and  demonstrations  are 
provided  using  Mountaintop  data  where  appropriate.  We  discuss  how  to  incorporate  Doppler  modulation 
generated  by  a  moving  receiver  as  well  as  that  generated  by  the  moving  emitter.  We  also  consider  the 
localization  of  multiple  emitters.  In  addition,  various  emitter  waveforms  are  examined  including  pulsed, 
periodic,  and  self-correlated  waveforms.  These  waveforms  often  cause  problems  for  localization  schemes 
because  multiple  peaks  arise  in  auto-  and  cross-correlation  functions.  The  multiple  peaks  may  generate 
false  emitter  locations  or  ghosts.  However,  we  shall  see  that  our  out-of-plane  multipath  technique  offers 
some  immunity  to  these  problems  because  it  uses  FDOA  (frequency  difference  of  arrival)  as  well  as 
TDOA  information  in  the  multipath  signals. 


12  MOVESfG  RECEIVER  AND  EMITTER 

If  the  receiver  is  also  moving,  the  Doppler  shifts  due  to  receiver  motion  must  be  considered  when 
determining  the  FDOA.  The  total  Doppler  shift  for  a  terrain-scattered  signal  will  be  the  sum  of  the 
Doppler  shift  due  to  each  moving  platform.  This  is  illustrated  in  Figure  7.1  which  shows  isodoppler 
contours  for  the  moving  receiver,  the  moving  emitter,  and  the  moving  receiver  and  emitter.  We  will  show 
that  the  addition  of  receiver  motion  to  the  moving  emitter  localization  problem  requires  only  a  simple 
Doppler  correction  when  the  receiver  velocity  and  heading  are  assumed  to  be  known. 

In  Appendix  C  the  differential  Doppler  was  shown  to  be 


f  =fs-fdp 


(7.1) 


where  R^is  the  unit  vector  along  R,  directed  toward  the  receiver,  and  R,is  the  unit  vector  along  R, 
directed  toward  the  scatterer.  The  terms  fs  and  fdp  denote  the  Doppler  shifts  for  the  scattered  and 
direct-path  signal,  respectively.  The  bistatic  triangle  is  shown  in  Figure  7.2. 

The  additional  Doppler  shift  caused  by  the  receiver  motion  can  be  included  in  (7.1)  as 

f  f  s  f  dp  ^  2^ 

=  (v,.r,-v,r„)+Kr,-v,.r„} 

where  the  term  in  parentheses  is  the  original  expression  and  the  term  in  brackets  is  the  contribution  due 
to  receiver  motion.  The  latter  term  is  known  completely  since  the  receiver  heading  and  velocity  are 


93 


Figure  7.1.  Doppler  contours  for  (a)  moving  receiver,  (b)  moving  emitter,  and  (c)  moving  receiver  and  emitter 


Scatterer 


Figure  7.2.  Bistatic  triangle  formed  by  the  receiver,  emitter,  and  scatterer. 
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assumed  to  be  known,  and  the  scatterer  and  emitter  azimuthal  angles  are  measured.  Hence,  the  results  of 
Chapters  3-5  can  be  directly  extended  to  the  moving  emitter  case.  We  note  that  for  a  specific  receive 
beam,  this  term  will  be  constant  since  the  Doppler  shift  due  to  the  receiver  is  given  hyVjX  cos  ^ . 


7.3  PULSED  EMITTER  WAVEFORM 

In  this  subsection  we  localize  a  periodic,  pulsed  emitter  using  the  estimator  designed  for 
inhomogeneous  clutter.  Up  to  this  point,  the  emitter  waveforms  were  assumed  to  be  white.  The  pulsed 
emitter  was  an  airborne  radar  transmitter  used  to  collect  bistatic  clutter  data  for  reflectivity 
measurements.  Thus,  localizing  this  emitter  is  essentially  equivalent  to  localizing  passively  an  airborne 
radar  system.  Data  set  HOT-6138  will  be  used,  and  the  location  of  this  emitter  was  shown  in  Figure  6.2. 

The  transmitter  waveform  is  pulsed  CW  with  a  pulse  width  of  13  ps,  a  PRF  of  300  Hz,  and  a  peak 
power  of  2  kW.  The  duty  cycle  of  this  waveform  is  0.39  percent  and  the  average  power  is  only  7.8  Watts. 
The  signal  in  the  direct-path  beam  and  two  multipath  receive  beams  is  shown  in  Figure  7.3.  The  strongest 
signal  present  in  the  multipath  beams  is  the  emitter  signal  that  is  present  in  the  sidelobes;  however, 
scattered  energy  is  evident  out  to  about  200  ps. 

The  scatterer  parameters  of  delay,  Doppler,  azimuth,  and  SNR  were  estimated  using  the  same 
structure  that  was  used  for  the  CW  waveforms.  That  structure  was  shown  in  Figure  6.13.  The  number  of 
scatterers  detected  was  40,  and  the  location  of  the  scatterers  is  shown  in  Figure  7.4.  This  is  one  half  the 
number  detected  for  HOT-6132  which  has  nearly  the  same  geometry.  This  reduction  in  the  number  of 
scatterers  is  probably  due  to  the  reduction  in  average  power  and  loss  of  resolution  for  the  pulsed 
waveform.  The  localization  results  are  listed  in  Table  7.1.  The  error  in  the  range  estimates  is  16  percent 
which  is  much  larger  than  the  examples  in  Chapter  6.  If  the  tme  emitter  edtitude  of  8  km  is  supplied  to  the 
estimator  so  that  only  range,  heading,  and  velocity  are  estimated,  then  the  results  are  188.6  km,  143.2, 
and  207.5,  respectively.  The  range  error  reduces  from  16  percent  to  only  3  percent.  It  is  interesting  that 
such  a  small  change  in  altitude  (from  11.8  to  8  km)  produced  such  a  large  change  in  range  estimate. 
Some  possible  explanations  for  this  behavior  were  given  in  Section  6.3.3. 

TABLE  7.1 

Results  for  Data  Set  HOT-6138 


Emitter 

Parameter 

True 

Value 

Estimated 

Value 

Error 

Value  CRLB 

CRLB 

a 

Range  (km) 

183.3 

153.2 

30.1 

23g 

1.3 

Heading  (deg) 

143.5 

141.9 

1.6 

1.0a 

1.5 

Velocity  (m/s) 

191.5 

149.5 

42 

49a 

0.85 

Altitude  (km) 

8 

11.8 

3.8 

2.2a 

1.7 

95 


Figure  7.4.  Map  ofscatterers  used  for  the  localization  of  data  set  HOT-6138. 


7.4  SELF-CORRELATED  EMITTER  WAVEFORM 

We  now  consider  correlation  in  the  emitter  waveform.  An  emitter  waveform  that  is  correlated  with 
itself  can  generate  multiple  TDOA  peaks  during  cross-correlation  operations.  The  multiple  TDOA 
estimates  can  produce  false  targets  or  ghosts.  This  is  illustrated  graphically  in  Figure  7.5  where  three 
hyperbolas  are  seen  to  intersect  at  nine  points. 


SENSOR  1 


SENSOR  3 


Figure  7.5.  Correlation  in  the  emitter  waveform  can  produce  multiple  TDOA  estimates  which  can  generate  false 
emitter  positions  when  using  multisensor  localization  TDOA  techniques. 
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To  investigate  the  effect  of  self-correlation  for  the  techniques  described  in  this  report,  we 
introduce  self-correlation  to  data  set  HOT-6067  as  shown  in  Figure  7.6.  The  autocorrelation  function  for 
the  correlated  data  set  is  shown  in  Figure  7.7.  The  peaks  are  close  enough  to  generate  ghosts  at  realistic 
ranges. 


Figure  7.6.  Schematic  for  introducing  self  correlation. 


Figure  7.7.  Autocorrelation  function  of  self-correlated  emitter  waveform. 


98 


The  estimator  designed  for  homogeneous  clutter  was  used  for  this  data  set  and  the  results  are  shown 
in  Figure  7.8.  Part  (a)  shows  the  log-likelihood  function  for  the  original  data  set  and  part  (b)  shows  the  log- 
likelihood  function  for  the  correlated  data  set.  There  is  some  spreading  of  the  peak  but  no  false  localizations 
occur.  Ghosts  are  not  generated  because  the  Doppler  present  in  the  multipath  signal  corresponds  only  to  the 
true  emitta-  since  it  is  due  to  the  emitta’s  motion.  The  correlation  manifests  itself  as  increased  noise  which 
blurs  the  peak  of  the  likelihood  ftinction. 


JAMMER  RANGE  (km)  JAMMER  RANGE  (km) 

(a)  (b) 

Figure  7.8.  Log-likelihood  function  (a)  without  and  (b)  with  self-correlation.  No  false  targets  are  generated. 

7.5  MULTIPLE  EMITTERS 

The  localization  techniques  described  in  this  report  are  readily  extended  to  the  case  of  multiple 
emitta-s.  In  this  section  we  consider  three  possible  approaches  for  use  with  the  estimator  for  homogeneous 
clutter.  The  basic  scenario  is  shown  in  Figure  7.9.  This  figure  shows  a  single  multipath  beam  and  two 
reference  beams,  one  for  each  emitta:.  The  two  emitter  waveforms  are  assumed  to  be  uncorrelated. 


Emitter  1 


Reference  Beam  1 


Si(n)^ 

Receiver 

z(n)  ^ 

S 

S2(n) 

- -~< 

Multipath  Beam  1 


Reference  Beam  2 


Emitter  2 


Figure  7.9  Localization  of  two  emitters.  The  multipath  beam  contains  energy  from  both  emitters.  Separate  reference 
beams  collect  the  two  emitter  waveforms. 


7.5.1  Sub-optimum  Technique 

The  simplest  approach  is  to  localize  the  emitters  one  at  a  time  while  essentially  ignoring  the 
presence  of  the  other  emitters.  This  is  not  optimum  because  the  other  emitters  will  produce  TSI  in  the 
receive  beam  which  effectively  raises  the  noise  level  of  the  receiver.  The  situation  can  be  improved 
somewhat  by  assuring  that  the  receive  and  reference  beams  are  orthogonal.  This  assures  that  unwanted 
emitter  energy  does  not  enter  through  the  sidelobes.  An  example  using  this  sub-optimum  technique  with 
Mountaintop  data  will  be  given  in  Section  7.5.4. 


7.5.2  Full  Maximum  Likelihood  Approach 

The  maximum  likelihood  solution  described  in  Chapter  4  can  be  extended  to  multiple  emitters. 
The  cost  is  fairly  high,  however.  For  example,  to  localize  two  emitters,  the  dimensions  of  the  parameter 
search  doubles,  as  does  the  size  of  the  covariance  matrix  that  is  inverted  at  each  step  of  the  search. 

To  extend  the  ML  solution  to  two  emitters,  Si  and  S2,  we  rewrite  (3.21)  as 
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and  solve  a  =  (s^s)  "s^zfor  the  filter  coefficients.  The  length  of  the  coefficient  vector  a  will  be  IP, 
where  P  is  the  length  of  the  original  filter.  The  parameter  vector  is  written  as 

e  =  {R„H„V„A„R2,H2,V^,A^}.  .  (7.4) 

To  maximize  the  likelihood  function  it  is  necessary  to  search  over  the  entire  parameter  vector  for  the 
minimum  output  of  the  loc2ilization  filter. 

7.5.3  Maximum  Likelihood  EM  Approach 

The  expectation-maximization  (EM)  algorithm  [74]  can  be  used  to  simplify  the  full  ML  approach 
described  above.  We  write  the  received  signal  for  two  emitters  with  parameter  vectors  6\  and  &i  as 

z{n)  =  mQ^{n)  +  mg^{n)  +  g{n)  ,  (7.5) 

where  g(n)  is  additive  Gaussian  noise.  Using  the  current  estimate  for  one  of  the  emitters,  di  say,  we  obtain 
Zj  (ji)  =  z{n)  -  m  -  in)  .  (7.6) 

This  is  an  estimate  of  the  signal  due  to  emitter  2  alone  given  that  m-  (n)  is  the  signal  due  to  emitter  1. 

We  then  use  the  estimated  signal  Z2in)  to  estimate  ft  as  we  did  in  Chapter  3  when  only  one  emitter  was 
present.  Then  the  estimate  of  the  signal  due  to  emitter  1  is  modified  to  be 

Ij  (n)  =  z(n)  -  m (n)  (7.7) 

and  proceed  to  refine  our  estimate  of  ft.  We  continue  this  iterative  procedure  until  the  parameter  sets 
converge.  Initial  guesses  at  the  parameter  sets  can  be  obtained  using  the  sub-optimal  method  described 
above. 
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7.5.4  Two  Emitter  Example  Using  the  Sub-optimum  Technique 

In  this  section  we  will  localize  two  emitters  in  Mountaintop  data  using  the  sub-optimum  approach 
described  in  Section  7.5.1  The  two-emitter  data  set  will  be  generated  by  combining  two  separate  data 
sets,  each  of  which  contains  a  single  emitter.  This  approach  is  necessary  because  the  Mountaintop 
experiments  did  not  use  more  than  one  airborne  emitter.  The  Milligan  Gulch  data  sets  HOT-6161  and 
HOT-6165  were  used.  HOT-6161  was  localized  in  Chapter  6  and  the  geometry  was  shown  in  Figure  6.2. 
The  data  for  HOT-6165  was  collected  a  few  minutes  after  HOT-6161  while  the  emitter  was  traveling 
along  approximately  the  same  track. 

An  advantage  of  studying  multiple  emitters  with  a  data  set  constructed  in  this  manner  is  that  one 
can  first  investigate  the  localization  of  the  emitters  separately,  and  then  investigate  the  performance  for 
the  combined  data  set.  The  conventional  and  MVDR  spatial  power  spectra  are  shown  in  Figure  7.10  (a) 
and  (b)  for  the  individual  data  sets  HOT-6161  and  HOT-6165,  respectively.  The  spatial  spectra  for  the 
combined  data  set  are  shown  in  (c).  The  conventional  spatial  power  spectral  estimate  does  not  resolve  the 
two  emitters  but  the  MVDR  spectrum  does  due  to  its  high  resolution  properties.  The  eigenvalues  of  the 
spatial  covariance  matrix  are  shown  in  part  (d).  The  two  strong  eigenvalues  correspond  to  the  strong 
direct  paths  and  the  weaker  eigenvalues  correspond  to  the  emitter  multipath.  Note  that  even  the  weakest 
eigenvalue  is  a  full  10  dB  above  receiver  noise. 

The  sub-optimum  approach  localizes  an  emitter  in  the  presence  of  TSI  from  other  emitters.  The 
TSI  from  other  emitters  acts  to  increase  the  noise  level  and  will  reduce  localization  performance.  Some 
attempts  can  be  made  to  minimize  the  effect  of  the  other  emitters,  however.  For  example,  spatially 
adaptive  antenna  beams  can  be  used  for  the  reference  beams.  This  assures  that  emitter  reference  signal 
for  one  emitter  does  not  contain  any  energy  from  the  other  emitter.  The  adapted  antenna  patterns  for  the 
two  emitters  are  shown  in  Figure  7.11.  Sidelobe  cancellers  are  used  to  remove  from  the  multipath  beam 
the  direct-path  and  short-time-delay  signals  from  each  of  the  two  emitters. 

The  emitters  were  localized  and  the  results  are  shown  in  Figure  7.12.  Parts  (a)  and  (b)  are  the 
results  of  the  single  emitter  localization,  and  parts  (c)  and  (d)  are  for  the  two-emitter  data  set.  Both 
emitters  have  been  localized,  but  the  dynamic  range  of  the  log-likelihood  function  has  been  reduced 
significantly  as  the  result  of  the  apparent  noise  increase  due  to  the  presence  of  the  second  emitter. 
Nevertheless,  the  localization  is  successful  and  could  be  used  to  initiate  an  iterative  procedure  such  as  the 
EM  algorithm. 
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Figure  7 AO.  Conventional  and  MVDR  spatial  power  spectra  in  dB  above  thermal  noise  for  single-emitter  data  sets 
(a)  HOT-616  f  (b)  HOT-6165,  and  (c)  the  two-emitter  combination  of  the  data  sets.  The  dashed  lines  represent  the 
true  emitter  angles.  Part  (d)  shows  the  eigenvalues  for  the  combined  data  set. 


(a) 


(b) 


Figure  7.11.  Adapted  beam  patterns  for  the  two  reference  beams.  Spatial  nulls  are  placed  at  the  location  of  the 
other  emitter. 


7.6  CHAPTER  SUMMARY 

In  this  chapter  we  have  outlined  methods  by  which  the  basic  results  of  Chapters  3-6  can  be 
extended  to  include  a  broader  range  of  scenarios  and  class  of  emitter  signals.  We  showed  that  receiver 
motion  could  be  handled  by  a  simple  adjustment  to  the  Doppler  terms  for  the  moving  emitter.  We  also 
considered  the  pulsed  and  self-correlated  waveforms,  as  well  as  multiple  emitters,  and  demonstrated 
good  performance  using  Mountaintop  data  for  each  case.  In  the  next  chapter  we  present  conclusions  as 
well  as  some  additional  topics  for  further  research. 
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Figure  7.12.  Single  and  multiple  emitter  localization  results  for  data  set  HOT-6161  and  HOT-6165.  Parts  (a)  and 
(b)  show  the  single  emitter  results,  and  (c)  and  (d)  show  the  multiple  emitter  results  for  HOT-6161  and  HOT-6165, 
respectively. 
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8.  SUMMARY  AND  CONCLUSIONS 


8.1  SUMMARY 

The  work  described  in  this  report  established  how  to  localize  passively  a  moving  emitter  using  out- 
of-plane  multipath  signals  reflected  by  the  terrain.  This  is  a  novel  localization  technique  that  assumes  no 
a  priori  knowledge  on  the  location  of  the  multipath  sources.  The  emitter  parameters  of  range,  heading, 
velocity,  and  altitude  are  estimated  by  exploiting  the  correlation  between  the  direct-path  signal  and  the 
delayed  and  Doppler  modulated  multipath  signals.  In  this  section  we  summarize  the  results  and  suggest 
areas  for  future  research. 

Two  basic  assumptions  about  the  scattering  properties  of  the  terrain  lead  to  two  different 
maximum  likelihood  (ML)  estimators  of  emitter  parameters.  The  first  assumption  is  that  the  terrain 
scattering  is  fundamentally  homogeneous;  in  this  case  the  ML  estimator  is  found  to  have  the  structure  of 
a  time-varying  FIR  filter.  The  time-varying  nature  of  the  filter  is  due  to  the  Doppler  modulators  that  are 
computed  deterministically  for  each  trial  parameter  set.  The  trial  parameters  are  selected  from  a  dense 
multidimensional  parameter  grid  over  which  we  search  for  the  true  emitter  parameters.  The  filter 
coefficients  (which  are  nuisance  parameters)  are  computed  using  the  method  of  least  squares  at  each 
point  in  the  search.  The  trial  parameter  set  that  minimizes  the  filter  output  is  the  ML  estimate. 

The  second  assumption  is  that  the  terrain  scattering  is  fundamentally  inhomogeneous  and 
dominated  by  a  number  of  discrete  scatterers.  This  assumption  leads  to  a  two-part  estimator  which  first 
estimates  the  scattering  parameters  of  azimuth,  differential  delay,  and  Doppler,  and  then  estimates  the 
emitter  parameters  of  range,  heading,  velocity,  and  altitude.  The  nonlinear  localization  equations  are 
linearized  by  using  a  truncated  Taylor  series  expansion  and  the  emitter  parameters  are  estimated  by 
performing  a  coarse  least-squares  search  followed  by  an  iterate  minimization  technique. 

The  Cramer-Rao  lower  bounds  for  each  estimator  are  derived  and  used  to  study  estimator 
performance  for  a  number  of  scenarios.  The  effect  of  varying  the  emitter  and  several  receiver  parameters 
was  considered.  We  saw  that  estimator  performance  was  related  to  emitter  velocity,  as  well  as  other 
factors  such  as  emitter  power,  heading,  altitude,  etc.  Monte-Carlo  simulations  are  performed  to  validate 
the  bounds  calculations  and  overall  implementation  of  the  estimators  prior  to  their  use  with  field  data. 
We  also  discussed  potential  limitations  that  could  arise  when  applying  these  bounds  to  real  data.  These 
limitations  resulted  from  simplifications  in  the  models  that  could  be  accounted  for  by  the  introduction  of 
model  noise  above  the  thermal  noise  level.  The  model  noise  was  found  to  be  20  dB  for  data  set 
HOT-6067. 

In  Chapter  6  we  demonstrated  the  estimators  using  field  data  collected  at  White  Sands  Missile 
Range  during  the  DARPA/Navy  Mountaintop  program.  Several  emitters  were  localized  using  both 
estimators.  In  general,  the  estimator  designed  for  inhomogeneous  clutter  performed  better  than  the 
estimator  for  homogeneous  clutter.  This  difference  in  performance  can  probably  be  attributed  to  the 
inhomogeneous  nature  of  the  bistatic  clutter  at  WSMR.  However,  the  estimator  for  homogeneous  clutter 
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performed  better  than  the  estimator  for  inhomogeneous  clutter  when  the  bandwidth  of  the  receiver  was 
reduced  to  20  kHz  or  less.  For  such  low  bandwidths,  the  discrete  scattering  model  breaks  down  due  to  the 
loss  of  range  resolution. 

In  Chapter  7  we  showed  how  the  basic  results  could  be  extended  to  include  a  broader  range  of 
scenarios  and  class  of  emitter  signals.  We  discussed  the  incorporation  of  Doppler  modulation  generated 
by  a  moving  receiver  as  well  as  the  moving  emitter.  This  was  found  to  require  only  a  minor  change  in  the 
way  relative  Doppler  shift  is  computed.  We  also  considered  the  localization  of  multiple  emitters  and  used 
a  sub-optimum  technique  to  localize  two  enutters  simultaneously.  In  addition,  we  considered  various 
emitter  waveforms  including  pulsed  and  self-correlated  waveforms  and  provided  an  example  for  each 
using  Mountaintop  data.  In  both  cases  the  emitters  were  localized  without  any  additional  false  targets  (or 
ghosts). 

8.2  FUTURE  RESEARCH 

There  are  several  potential  areas  for  further  research  related  to  the  topics  in  this  report.  They  can 
be  roughly  divided  into  two  areas:  further  investigations  into  the  bistatic  scattering  properties  of  the 
terrain  and  further  work  on  the  estimation  techniques.  We  discuss  a  few  examples  in  this  section. 

A  very  important  question  is  whether  or  not  this  localization  technique  can  be  used  over  the  ocean. 
Sea  clutter  presents  two  potential  problems.  The  first  is  that  sea  clutter  exhibits  intrinsic  motion  that  will 
generate  additional  Doppler  shifts.  This  may  reduce  localization  performance  or  bias  estimates  if  not 
accounted  for.  The  second  problem  is  that  sea  clutter  may  be  too  homogeneous  to  use  the  estimator 
designed  for  inhomogeneous  clutter  unless  sea  spikes  are  prevalent.  In  addition,  the  general  bistatic 
scattering  levels  for  sea  clutter  may  be  much  lower  than  for  land  clutter;  however,  we  saw  in  Section  6.4 
that  there  is  considerable  excess  signal-to-noise  ratio  available  for  the  Mountaintop  parameters  and, 
based  on  that  observation,  we  believe  the  estimator  designed  for  homogeneous  clutter  will  work  over  the 
ocean.  This  assertion  should  be  validated  with  field  data,  however. 

There  are  several  issues  related  to  the  estimator  for  the  inhomogeneous  scattering  model  that 
should  be  investigated  further.  The  work  here  assumed  that  all  the  scatterers  that  were  detected  had 
variances  inversely  proportional  to  their  SNR.  However,  we  saw  in  Section  6.3.3  that  by  excluding  one  of 
eighty  scatterers,  the  localization  estimates  improved  significantly.  We  believe  the  parameters  for  this 
scatterer  have  been  corrupted  by  coherent  interference.  Future  work  could  investigate  effective  ways  of 
selecting  a  minimum  subset  of  scatterers  that  will  produce  good  localization  estimates.  An  alternative  to 
excluding  particular  scatterers  could  be  to  treat  them  as  nuisance  parameters  to  be  estimated  along  with 
the  emitter  parameters.  A  related  issue  is  the  effect  of  “false  scatterers”  on  the  localization  estimator.  For 
example,  a  noise  spike  or  bistatic  scattering  from  a  target  such  as  a  moving  aircraft  will  have  erroneous 
scattering  parameters  when  viewed  as  bistatic  clutter,  and  these  false  scatterers  could  degrade 
localization  estimates.  Future  work  should  include  a  study  of  the  behavior  of  this  estimator  in  the 
presence  of  erroneous  scattering  parameters.  The  Hough  transform  described  in  Chapter  2  may  be  a 
useful  tool  for  evaluating  graphically  the  contributions  of  individual  scatterers. 
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The  Doppler  terms  computed  in  this  report  were  based  on  a  flat-earth  approximation.  Future  work 
could  include  a  reformulation  of  the  geometry  using  site-specific  digital  elevation  data.  This  improved 
geometric  model  may  reduce  localization  errors  or  biases,  especially  for  mountainous  terrain  such  as  that 
found  at  WSMR. 

The  use  of  a  priori  scattering  information  was  not  considered  in  this  report;  however,  if  a  priori 
knowledge  of  the  location  of  strong  scatterers  was  available,  it  could  be  used  to  localize  emitters.  This 
may  permit  the  localization  of  stationary  emitters  that  would  otherwise  not  be  possible  using  the 
techniques  described  in  Chapters  3  and  4.  Sources  of  a  priori  knowledge  on  bistatic  scattering  could 
include  site-specific  bistatic  clutter  models  that  use  digital  elevation  data.  Another  source  of  terrain 
knowledge  could  be  extracted  from  the  monostatic  clutter  data  collected  in  an  active  mode.  The  idea 
would  be  to  compare  azimuth-delay  maps  measured  using  cross-correlation  of  the  emitter  waveform  to 
maps  generated  based  on  trial  emitter  locations  and  the  a  priori  clutter  information. 

We  close  this  section  by  noting  that  in  this  report  we  have  demonstrated  the  ability  to  localize 
moving  emitters,  using  out-of-plane  multipath,  for  one  system  operating  at  one  site.  The  Mountaintop 
radar  has  a  UHF  phased-array  receiver  and  was  operating  over  desert  and  mountainous  terrain  at  WSMR. 
To  fully  explore  the  usefulness  of  these  localization  techniques  for  other  applications,  modeling  and 
measurement  programs  should  be  carried  out  for  additional  systems  and  over  additional  terrain  types. 
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APPENDICES 


Appendices  A-D  contain  the  geometrical  relations  and  derivations  necessary  to  support  the  work 
contained  in  this  report.  Included  is  the  determination  of  the  Doppler  shift  terms  for  two-  and 
three-dimensional  approximations  as  well  as  the  general  three-dimensional  solution.  Also  included  are 
the  differential  terms  necessary  to  form  the  Fisher  information  matrix. 
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APPENDIX  A 

DOPPLER  TERMS  FOR  TWO-DIMENSIONAL  APPROXIMATION 


The  two-dimensional  approximation  to  the  bistatic  geometry  of  concern  here  makes  use  of  the  fact 
that  the  platform  ranges  are  generally  much  greater  than  the  platform  altitudes.  For  example,  the  height  of 
the  Mountaintop  receiver  above  the  ground  is  approximately  1.2  km  and  the  emitter  altitude  is  7-8  km 
while  the  range  to  the  emitter  is  typically  100-200  km.  Therefore,  we  shall  assume  for  this  section  that 
the  receiver,  emitter,  and  scatterer  are  all  located  in  the  same  horizontal  plane.  The  top-down  view  of  the 
bistatic  triangle  is  shown  in  Figure  A.l.  We  have  also  assumed  that  the  earth  is  flat  and  defined  the 
azimuthal  angle  of  the  emitter  relative  to  the  receiver  to  be  zero. 


Scatterer 


Figure  A.l.  Bistatic  triangle  formed  by  the  receiver,  emitter,  and  scatterer. 


The  goal  is  to  obtain  an  expression  for  Doppler  shift  /as  a  function  of  the  emitter  parameters  (R, 
H,  V),  and  the  measurable  quantities  of  bistatic  delay  D  and  scatterer  azimuth  0.  The  bistatic  delay  D  is 
related  to  the  range  to  the  emitter  R,  the  range  to  the  scatterer  /?„  and  the  range  firom  the  emitter  to  the 
scatterer  /?,  as 

D  =  R^+R,-R.  (A.1) 

Using  the  law  of  cosines  we  can  solve  for  the  range  to  the  scatterer  Rr  in  terms  of  the  emitter  range  R, 
bistatic  delay  D,  and  angle  0 

R  ^_R£±dV2_  (^2) 

D  +  R-Rcos^ 

The  range  from  the  emitter  to  the  scatterer  /?,  is  obtained  by  rewriting  (A.l)  as 
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R,  =R  +  D-R^  . 

The  angle  yean  be  obtained  using  the  law  of  sines. 


y  =  sin' 


R^  sin  <!> 


and  the  angle  yf  can  be  written  in  terms  of  H  and  y, 

y/  =  180-H-y  . 

The  Doppler  shift  observed  at  the  location  of  the  scatterer  due  to  the  emitter  motion  is 


fs  =jCOS\{/  , 

which  can  also  be  written  as 


fs  =~cos{H  +  y)  , 

where  we  have  made  use  of  (A.5).  The  Doppler  shift  for  the  direct-path  signal  is 

U=jCOs{lS0-H) 


The  stationary  receiver  measures  the  relative  Doppler  shift  between  the  direct  and  scattered  paths, 

f  =  fs-fdp 

=  y(-  cos{H  +  y)+  cos(/f)) . 

A 

The  angle  yean  be  expressed  in  terms  of  the  desired  quantities  by  using  (A.2)  and  (A.3)  in  (A.4)  as 
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7 


=  sin 


sin  (j) 
~~R. 


/ 


=  sin 


-1 


^  RD  +  D^/l  V  ^  ^ 
'  Ism  (p 


D  +  R-Rcos<j) 


f 


R^D- 


D  +  /?-/?cos^ 


=  sin 


-1 


_ (d^ +2RD)sm(f>  ^ 

+  2/?^  +  2RD  -  2(r^  +  DR)cos<I> 


and  finally. 


r 


cos{h)  -cos 


r 


V 


H  +  sin 


-1 


( 

+2/2Z))sin^ 

DI  +  2 

\[R^+RDp 

\ 

o 

o 

which  is  the  final  result. 


(A.  10) 


(A.  11) 


To  summarize  this  appendix,  we  have  derived  an  expression  for  the  differential  Doppler  shift, 
between  the  scattered  and  direct  path,  in  terms  of  the  emitter  parameters  (R,  H,  V)  and  the  measured 
parameters  (D,  0).  This  equation  was  derived  based  on  the  assumption  that  the  receiver,  scatterer,  emitter, 
and  emitter  velocity  are  contained  in  the  same  plane.  Appendix  B  provides  an  approximation  for  the  case 
when  the  emitter  altitude  is  nonzero  and  Appendix  C  treats  the  general  case. 
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APPENDIX  B 

DOPPLER  TERMS  FOR  THREE-DIMENSIONAL  APPROXIMATION 


The  development  in  Appendix  A  contains  no  provision  for  emitter  altitude  estimation  because  it 
assumes  that  both  platforms  are  at  zero  altitude.  In  this  appendix,  we  relax  this  restriction  slightly  to 
allow  the  emitter  to  have  a  nonzero  altitude  that  is  small  compared  to  the  range  terms.  The  Doppler  shift 
at  the  scatterer  for  the  two-dimensional  case  was 


f^=jCOSY  ,  (B-1) 

where  yrwas  the  azimuth  angle  measured  from  the  emitter  velocity  vector  to  the  scatterer.  If  the  emitter  is 
elevated,  then  Xj/  will  be  a  conical  angle  measured  relative  to  the  emitter  velocity  vector  to  the  scatterer. 
Using  direction  cosines,  we  can  write 

cos  =  cos  cos  ,  (®-2) 

where  and  £„  are,  respectively,  the  azimuth  and  elevation  angles  to  the  scatterer.  We  can  write  the 
cosine  of  the  elevation  angle  in  terms  of  the  scatterer  range  R,  and  emitter  height  h,  as 


cosEj  = 


(B.3) 


and  rewrite  (A.9)  as 


f 

(/  7  2  \  ^ 

11 

_ 

1 

‘  n 

cos(H  +  7)-!- cos  (H) 


and  repeat  the  algebra  in  (A.  10)  and  (A.1 1).  The  final  equation  is 


/  = 


V 


cos 


(H)- 


R 


cos 


sin  ' 


{d^ +2RD)sin<t) 

dI+2{r^+RdII\-cos(I))^ 


(BA) 


(B.5) 


with  R,  given  by 


^  +3D^/2  +  DR-{r^ +Dr)cos(I> 

'  D  +  R-Rcos(l> 

This  assumes  the  height  of  the  emitter  is  not  large  enough  to  change  significantly  the  range  terms  R  or  R,. 
Appendix  C  treats  the  general  case  where  no  such  assumptions  or  restrictions  are  necessary. 
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APPENDIX  C 


DOPPLER  TERMS  FOR  COMPLETE  THREE-DIMENSIONAL 

SOLUTION 

In  this  appendix  we  consider  the  general  three-dimensional  case  shown  in  Figure  C-1.  The  receiver 
is  located  at  the  origin  of  the  rectangular  coordinate  system,  and  the  conical  angles  to  the  scatterer  C„ 
and  emitter  are  measured  relative  to  the  jc-axis.  We  desire  to  obtain  an  expression  relating  the  Doppler 
shift /to  the  emitter  parameters  (R,  H,  V,  ht),  and  the  measurable  quantities  of  bistatic  delay  D,  scatterer 
cone  angle  C„,  and  emitter  cone  angle 


^  t  -Ers  s 

✓  / 


Figure  C.l.  General  three-dimensional  bistatic  geometry. 


An  expression  for  the  receiver  angle  between  the  emitter  and  the  scatterer,  can  be  obtained 
using  direction  cosines  [47] 

cos 7^  =  cos £„  cos E„  cos(>l„  -  )  +  sin  sin  ,  (C.l) 

where  the  angles  are  as  shown  in  the  figure.  If  7  were  known,  we  could  proceed  as  we  did  above  in 
(A.2)-(A.ll);  however,  this  is  not  the  case  so  we  shall  proceed  as  follows.  We  first  write  the  measured 
conical  angles  in  terms  of  azimuth  and  elevation  angles, 

cos7„  =cosA„cos£'^,  (C.2) 

and 

cos7„  =cosA„cos£'„  .  (C.3) 
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The  receiver  azimuth  and  elevation  angles  to  the  emitter  can  be  written  as 


.  _if  cosC„ 

A„  =  sin  - - 

cosE, 


E„  =  sin- 


if 


Similarly,  the  receiver  azimuth  and  elevation  angles  to  the  scatterer  can  be  written  as 


A„  =  sin' 


=  sin" 


cos  C^. 
cos  E_ 


[K-hl 


We  can  now  insert  (C.l)  into  (A.2)  to  yield 

„  RD-D^ll 

R  = - i _  g\ 

D  +  R  —  R(cos  E^^  cos  E^  cos(>i^^  “  )  +  sin  E^^  sin  E^  ) 

and  solve  numerically  using  Newton’s  method.  The  numerical  method  is  required  because  Rr  appears  on 
both  sides  of  (C.8)  because  of  the  elevation  angle  described  by  (C.7).  For  the  initial  values  required  to 
start  the  iterative  technique,  we  set  the  two  elevation  angles  to  zero,  and  then  (C.8)  reduces  to  the  two- 
dimensional  case  (A.2)  which  can  be  solved  directly. 

We  proceed  in  a  similar  manner  to  obtain  the  azimuth  and  elevation  angles  of  the  scatterer  relative 
to  the  emitter.  The  range  from  the  emitter  to  the  scatterer  was  given  by  (A.3)  as 

R,=R  +  D-R^  .  (A.3) 

The  elevation  eingle  to  the  scatterer  is 


E,  =  sin  ^ 

R. 


To  find  the  azimuth  angle  we  rewrite  (C.l)  for  the  transmitter  as 


cos  Y,  =  cos  cos  cos(A,  -  )  -t-  sin  sin  E^ 


(C.IO) 


and  solve  for  A,.  The  angle  y,  is  determined  using  the  law  of  cosines 


(C.ll) 


r, 


cos 


^R^+Rf-R^  ^ 


With  the  bistatic  triangle  determined,  we  now  proceed  to  determine  relative  Doppler  shift.  For  this 
development  it  will  be  convenient  to  use  vectors  and  dot  products.  As  a  final  check  on  our  derivations, 
we  shall  show  that  the  three-dimensional  case  will  reduce  to  the  two-dimensional  case  described  above 
when  the  receiver  and  emitter  altitudes  are  set  to  zero. 


We  now  introduce  two  unit  vectors:  is  the  unit  vector  along  R,  directed  toward  the  transmitter 

and  R,  is  the  unit  vector  along  R,  directed  toward  the  scatterer.  Expressions  for  these  vectors  are 


R^j  =  X  cos  cos  E^+y  sin  A^  cos  +  z  sin 

(C.12) 

R,  =  xcosA,  cosE, -l-ysin  A,  cosE, -t-zsinE,  . 

The  emitter  velocity  vector  is 


V  =  V(xcos/f +  ysin/f). 
The  differential  Doppler  shift  is 


(C.13) 


J  J  s  J  dp 

=  V  R,-V  R, 

=  V.(r,+R„) 

where  we  have  made  use  of  the  fact  that  R^  R,r  •  After  plugging  in  terms  we  can  write 


(C.14) 


/  = — (xcosH-hysin/f)-[x(cosA,  cosE,  -l-cosA„  cosE„)-l- 

A 

y(sin  A,  cos  E,  -I-  sin  A„  cos  E„  )-l-  z(sin  E,  -  sin  E„  )]  . 


After  carrying  out  the  dot  product  we  obtain 


(C.15) 


(C.16) 


which  is  the  desired  expression.  All  the  terms  on  the  right  side  of  the  equation  depend  on  measured  or 
known  quantities  and  emitter  parameters  we  want  to  estimate.  (To  obtain  to  from/  we  multiply  by  2n.) 
As  a  check  on  this  equation,  we  set  the  receiver  height,  emitter  height,  and  Art  equal  to  zero  and  compare 
to  (A.9),  the  two-dimensional  approximation.  For  this  case,  (C.16)  yields 
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(C.17) 


/  =  ^{cosH{cosA,  +l)+sin//sin  a,) 

A 

V  t  \ 

=  —  ^cos  H  cos  A,  +  sin  H  sin  A,  +  cos  H  ) 

A 

=  Y {cos{H  -  A, )+  cos //) . 

A 

Since  A,,  the  angle  to  the  scatterer  from  the  emitter,  is  I8O-7,  where  7  was  defined  in  Figure  A.l,  we  can 
rewrite  this  result  as 


/  =  “180  +  7)+cos//) 

A 


=  Y  (~  cos(H  +  7)+ cos ) 
A 


(C.18) 


which  is  equivalent  to  (A.9). 

Delay-Doppler  ridges  were  computed  using  the  equations  derived  in  Appendices  A-C  and  the 
results  are  shown  in  Figure  C.2.  The  two-dimensional  treatment  used  to  derive  Equation  (A.ll)  over- 


Figure  C.2.  Delay-Doppler  ridge  for  three  methods  of  calculation.  The  emitter  range  is  100  km  and  its  heading  is 
90  degrees.  The  altitudes  for  the  emitter  and  receiver  are  8  km  and  1.2  km,  respectively. 
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predicts  the  Doppler  frequency  slightly  because  it  uses  azimuthal  angles  instead  of  conical  angles.  The 
approximate  three-dimensional  solution  that  resulted  in  Equation  (B.4)  is  very  close  to  the  full  three- 
dimensional  solution  because  it  accounts  for  the  emitter  conical  angles.  One  can  see  from  the  figure  that 
the  approximate  and  full  three-dimensional  solutions  are  almost  identical  for  this  case. 
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APPENDIX  D 

ELEMENTS  OF  THE  FISHER  INFORMATION  MATRIX 


In  this  appendix  we  derive  the  expressions  for  the  individual  elements  of  the  Fisher  information 
matrix  corresponding  to  the  statistical  formulation  that  was  presented  in  Chapter  4.  It  was  shown  there 
that  the  elements  of  the  Fisher  matrix  F  are 

dm^  (n)  dm^  (n)  ^  dm.  (n)  dm.  {^)  ^  ^  ^ 

dd,  de,  de,  ddj  ’ 

I  J  i  •'  / 

with 

m(n)-^apS(n-Dp)e^ 

p=i 

The  scattering  coefficient  Op  and  the  emitter  waveform  s(n)  are  complex.  The  real  and  imaginary  terms 
can  be  written  in  an  amplitude-phase  form  as 


P 

(«)  =  XK  ih"  “  ^1  -  -O/-  )  +  +  9a,  ) 

(D.3) 

p 

(n)  =  I  s(n  -  Dp  )|  sm{a)p  in-Dp)  +  (Ps(n-D,)  +9a,)  ■ 

/7=1 


For  convenience,  we  combine  the  phase  terms  in  (Pp  and  drop  the  absolute  value  symbols  by  redefining  ap 
and  s(n)  to  be  amplitudes  and  rewrite  (D.3)  as 


fn,(n)  =  '^aps(n-Dp  )cos(oDp  (0}i  +  (Pp) 

p=l 

m,. (n)  =  ^aps{n  - )sin(©p {e)n  +  (pp). 

p=\ 


(DA) 


The  dependence  of  O)  on  the  emitter  parameters  0has  been  indicated  by  their  inclusion  in  the  sum. 


The  entire  parameter  vector  contains  the  desired  emitter  parameters  and  the  2P  nuisance 
parameters.  This  is  written  as  {R,Hy,A, (p^  ,...,(pp ,a.^,...,ap } .  The  required  derivatives  are 

dm,(n)  dm.jn) 

<?6>,  ’  de. 
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dm^iri)  dm.(n) 
dm^(n)  dm.{n) 

^(Pp  ’  ^<Pp 

where  0 denotes  the  desired  emitter  parameters  {/?,  H,  V,  A}.  Differentiating  (D.4)  with  respect  to  the 
desired  parameters  yields 


^  D-\ 


To  complete  the  differentiation,  we  require  the  evaluation  of 

d(Op{e) 


(D.IO) 


for  each  of  the  desired  parameters.  For  range  R,  heading  H,  and  velocity  V,  we  differentiate  (A.  11)  with 
respect  to  each  parameter.  For  convenience,  we  first  rewrite  (A.l  1)  as 


VI  f  II 

f--;r  cos{h)-cos  H  +  sm~  — 

A  V 


(D.ll) 


/t  =  (Z)^  +  2RD)sm<l> 


(D.12) 


v  =  D^+2R^+  2RD  -2(r^+  Dr)cos  (p  . 


The  required  derivatives  are 


(D.13) 


doM 


dR  A  I  V 


v(2Dsin^)-/t(l-cos0X‘^^  +  2D) 


(D.14) 
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dm^(n) 

da. 


^s{n-D  )cos(fflp {0}i  +  (pp) 


(D.21) 


=s{n  - D  )sin(ft)p {0)n  +  (pp)  (0-22) 

da, 

for  the  amplitude  terms.  The  sum  over  P  is  not  present  as  it  was  for  (D.8)  and  (D.9)  because  we  have 
assumed  that  only  the  p‘^  term  in  the  summation  depends  on  a,  or 

With  the  above  derivatives  determined,  the  Fisher  information  matrix  is  formed  by  substituting  the 
above  expressions  in  (D.l).  The  results  so  far  contain  the  emitter  waveform  s(n)  explicitly.  If  we  assume 
that  the  emitter  waveform  has  variance  and  a  correlation  length  that  is  less  than  the  filter  coefficients 

spacing,  we  can  write  the  Fisher  matrix  in  terms  of  (7^ .  When  this  is  true,  the  expected  values  of 
products  of  the  form 
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(D.23) 


where  we  have  used  the  real-valued  notation  introdueed  for  convenience,  reduce  to 


(D.24) 


Hence,  for  the  emitter  parameters  terms,  the  diagonal  elements  of  the  Fisher  information  matrix  can  be 
written  as 


d(o^{e) 


and  the  general  expression  is 

•  cr"  ^  '  de,  d0j 

For  the  nuisance  parameters  on  the  diagonal,  the  phase  terms  reduce  to 


(D.25) 


(D.26) 


p  ^ 


and  the  amplitude  terms  reduce  to 


(D.27) 


(D.28) 


In  these  equations,  the  quantity  — |-is  the  SNR  for  the  direct-path  signal  and  the  quantity  - - — — 

represents  the  integrated  SNR  at  the  p**  filter  coefficient.  The  final  form  of  the  Fisher  matrix  for  this 
assumption  was  noted  in  Chapter  5  and  shown  to  contain  a  lot  of  zero  terms  off  the  diagonal. 
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